Artifact aa69303de31ed67aa3c040d6c94724e4636b1b8b1147adb76ddcb2206bf16031:
- Executable file
r37/packages/taylor/tayfns.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: 47511) [annotate] [blame] [check-ins using] [more...]
module TayFns; %***************************************************************** % % Simplification functions for special functions % %***************************************************************** exports taysimpexpt, taysimpatan, taysimplog, taysimpexp, taysimptan, taysimpsin, taysimpsinh, taysimpasin; imports % from the REDUCE kernel: !*f2q, !:minusp, addsq, aeval, denr, domainp, eqcar, evenp, freeof, invsq, kernp, lastpair, let, lprim, lnc, mk!*sq, mksq, multsq, mvar, negsq, neq, nlist, nth, numr, over, prepd, prepsq, quotsq, retimes, reval, reversip, simp, simp!*, simplogi, simplogsq, subs2!*, subsq, subtrsq, % from the header module: !*tay2q, !*TayExp2q, constant!-sq!-p, cst!-Taylor!*, find!-non!-zero, get!-degree, has!-TayVars, make!-cst!-powerlist, make!-Taylor!*, prune!-coefflist, set!-TayCoeffList, set!-TayFlags, set!-TayOrig, TayCfPl, TayCfSq, TayCoeffList, TayFlags, TayGetCoeff, Taylor!*p, Taylor!-kernel!-sq!-p, Taylor!:, TayMakeCoeff, TayOrig, TayTemplate, TayTpElNext, TayTpElOrder, TayTpElPoint, TayTpElVars, TayTpVars, TayVars, TpNextList, % from the module Tayintro: confusion, delete!-nth, delete!-nth!-nth, replace!-nth, replace!-nth!-nth, Taylor!-error, Taylor!-error!*, var!-is!-nth, % from the module Tayutils: addto!-all!-TayTpElORders, get!-cst!-coeff, is!-neg!-pl, smallest!-increment, subtr!-degrees, Taylor!*!-constantp, Taylor!*!-zerop, % from the module Taybasic: addtaylor, addtaylor!-as!-sq, invtaylor, makecoeffs0, makecoeffpairs, makecoeffpairs1, multtaylor, multtaylor!-as!-sq, multtaylorsq, negtaylor, negtaylor1, quottaylor, % from the module TayExpnd: taylorexpand, % from the module Taysimp: expttayrat, taysimpsq, taysimpsq!*, % from the module Taydiff: difftaylorwrttayvar, % from the module TayConv: prepTayCoeff, prepTaylor!*, % from the module Tayfrontend: taylorcombine, taylortostandard; fluid '(!*taylorkeeporiginal !*!*taylor!-epsilon!*!* frlis!*); symbolic procedure taysimpexpt u; % % Argument is of the form ('expt base exponent) % where both base and exponent (but a least one of them) % may contain Taylor kernels given as prefix forms. % Value is the equivalent Taylor kernel. % if not (car u eq 'expt) or cdddr u then confusion 'taysimpexpt else if cadr u eq 'e then taysimpexp {'exp, caddr u} else begin scalar bas, expn; bas := taysimpsq simp!* cadr u; expn := taysimpsq simp!* caddr u; if constant!-sq!-p bas then return taysimpexp {'exp,mk!*sq multsq(simp!*{'log,mk!*sq bas},expn)}; if null kernp bas then if not(denr bas = 1) then return mksq({'expt,prepsq bas,prepsq expn},1) else if domainp numr bas then return taysimpexp {'exp, prepsq multsq(simp!* {'log,prepd numr bas},expn)} else return mksq({'expt,prepsq bas,prepsq expn},1); if fixp numr expn and fixp denr expn then return !*tay2q expttayrat(mvar numr bas,expn); if denr expn = 1 and eqcar(numr expn,'!:rn!:) then return !*tay2q expttayrat(mvar numr bas,cdr numr expn); bas := mvar numr bas; return if Taylor!*p bas then if Taylor!-kernel!-sq!-p expn then if TayTemplate bas = TayTemplate mvar numr expn then taysimpexp {'exp, mk!*sq taysimpsq multtaylor!-as!-sq( expn, taysimplog {'log,bas})} else mksq({'expt,bas,mvar numr expn},1) else if not has!-TayVars(bas,expn) then begin scalar logterm; logterm := taysimplog{'log,bas}; return if Taylor!-kernel!-sq!-p logterm then taysimpexp{'exp, multtaylorsq(mvar numr logterm, expn)} else taysimpsq simp!* {'exp,mk!*sq multsq(logterm,expn)} end else mksq({'expt,bas,mk!*sq expn},1) else if Taylor!-kernel!-sq!-p expn then if not has!-TayVars(mvar numr expn,bas) then taysimpexp{'exp, multtaylorsq(mvar numr expn, simp!*{'log,bas})} else if Taylor!*!-constantp mvar numr expn then taylorexpand( simp!* {'expt,bas, prepTaylor!* mvar numr expn}, TayTemplate mvar numr expn) else mksq({'expt,bas,mk!*sq expn},1) else mksq({'expt,bas,mk!*sq expn},1); end; put('expt,'taylorsimpfn,'taysimpexpt); symbolic procedure TayCoeffList!-union u; if null cdr u then car u else TayCoeffList!-union2 (car u, TayCoeffList!-union cdr u)$ symbolic procedure TayCoeffList!-union2 (x, y); % % returns union of TayCoeffLists x and y % << for each w in y do if null (assoc (car w, x)) then x := w . x; x >>$ symbolic procedure inttaylorwrttayvar(tay,var); % % integrates Taylor kernel tay wrt variable var % inttaylorwrttayvar1(TayCoeffList tay,TayTemplate tay,var)$ symbolic procedure inttaylorwrttayvar1(tcl,tp,var); % % integrates Taylor kernel with TayCoeffList tcl and template tp % wrt variable var % Taylor!: begin scalar tt,u,w,singlist,csing; integer n,n1,m; u := var!-is!-nth(tp,var); n := car u; n1 := cdr u; tt := nth(tp,n); u := for each cc in tcl join << m := nth(nth(TayCfPl cc,n),n1); if TayTpElPoint nth(tp,n) eq 'infinity then << if m=1 then <<singlist := TayMakeCoeff( delete!-nth!-nth(TayCfPl cc,n,n1), TayCfSq cc) . singlist;nil>> else {TayMakeCoeff( replace!-nth!-nth(TayCfPl cc,n,n1,m-1), multsq(TayCfSq cc,invsq !*TayExp2q(-m + 1)))}>> else << if m=-1 then <<singlist := TayMakeCoeff( delete!-nth!-nth(TayCfPl cc,n,n1), TayCfSq cc) . singlist;nil>> else {TayMakeCoeff( replace!-nth!-nth(TayCfPl cc,n,n1,m+1), multsq(TayCfSq cc,invsq !*TayExp2q(m + 1)))}>>>>; w := {TayTpElVars tt,TayTpElPoint tt, if var member TayTpElVars tt then if TayTpElPoint tt eq 'infinity then TayTpElOrder tt - 1 else TayTpElOrder tt + 1 else TayTpElOrder tt, if var member TayTpElVars tt then if TayTpElPoint tt eq 'infinity then TayTpElNext tt - 1 else TayTpElNext tt + 1 else TayTpElOrder tt}; if singlist then begin scalar tpel; tpel := nth(tp,n); singlist := reversip singlist; if TayCfPl car singlist = '(nil) % no Taylor left then csing := TayCfSq car singlist else csing := !*tay2q make!-Taylor!*( singlist, replace!-nth( tp,n, {delete!-nth(TayTpElVars tpel,n1), TayTpElPoint tpel, TayTpElOrder tpel, TayTpElNext tpel}), nil,nil); csing := multsq(csing,simp!* {'log,nth(TayTpElVars tpel,n1)}) end; return (csing . make!-Taylor!*(u,replace!-nth(tp,n,w),nil,nil)) % % The following is not needed yet % % return make!-Taylor!*( % u, % replace!-nth(TayTemplate tay,n,w), % if !*taylorkeeporiginal and TayOrig tay % then simp {'int,mk!*sq TayOrig tay,var} % else nil, % TayFlags u) end; comment The inverse trigonometric and inverse hyperbolic functions of a Taylor kernel are calculated by first computing the derivative(s) with respect to the Taylor variable(s) and integrating the result. The derivatives can easily be calculated by the manipulation functions defined above. The method is best illustrated with an example. Let T(x) be a Taylor kernel depending on one variable x. To compute the inverse tangent T1(x) = atan(T(x)) we calculate the derivative T'(x) T1'(x) = ----------- . 2 1 + T(x) (If T and T1 depend on more than one variable replace the derivatives by gradients.) This is integrated again with the integration constant T1(x0) = atan(T(x0)) yielding the desired result. If there is more than variable we have to find the potential function T1(x1,...,xn) corresponding to the vector grad T1(x1,...,xn) which is always possible by construction. The prescriptions for the eight functions asin, acos, asec, acsc, asinh, acosh, asech, and acsch can be put together in one procedure since the expressions for their derivatives differ only in certain signs. The same remark applies to the four functions atan, acot, atanh, and acoth. These expressions are: d 1 -- asin x = ------------- , dx sqrt(1-x^2) d -1 -- acos x = ------------- , dx sqrt(1-x^2) d 1 -- asinh x = ------------- , dx sqrt(1+x^2) d 1 -- acosh x = ------------- , dx sqrt(x^2-1) d 1 -- atan x = --------- , dx 1 + x^2 d -1 -- acot x = --------- , dx 1 + x^2 d 1 -- atanh x = --------- , dx 1 - x^2 d 1 -- acoth x = --------- , dx 1 - x^2 together with the relations 1 asec x = acos - , x 1 acsc x = asin - , x 1 asech x = acosh - , x 1 acsch x = asinh - x . This method has one drawback: it is applicable only when T(x0) is a regular point of the function in question. E.g., if T(x0) = 0, then asech(T(x)) cannot be calculated by this method, as asech has a logarithmic singularity at 0. This means that the constant term of the series cannot be determined by computing asech(T(x0)). In that case, we use the following relations instead: asin z = -i log(i z + sqrt(1 - z^2)), acos z = -i log(z + sqrt(z^2 - 1)), 1 1 + i z atan z = ----- log ---------, 2 i 1 - i z -1 i z + 1 acot z = ----- log ---------, 2 i i z - 1 asinh z = log(z + sqrt(1 + z^2)), acosh z = log(z + sqrt(z^2 - 1)), 1 1 + z atanh z = --- log -------, 2 1 - z 1 z + 1 acoth z = --- log -------. 2 z - 1 These formulas are applied at the following points: infinity for all functions, +i/-i for atan and acot, +1/-1 for atanh and acoth. There are still some branch points, where the calculation is not always possible: +1/-1 for asin and acos, and consequently for asec and acsc, +i/-i for asinh, acosh, asech and acsch. In these cases, the above formulas are applied as well, but simplification of the square roots and the logarithm may lead to a rather complicated result; symbolic procedure taysimpasin u; if not (car u memq '(asin acos acsc asec asinh acosh acsch asech)) or cddr u then confusion 'taysimpasin else Taylor!: begin scalar l,l0,c0,v,tay0,tay,tay2,tp,singlist; tay0 := taysimpsq simp!* cadr u; if not Taylor!-kernel!-sq!-p tay0 then return mksq({car u,mk!*sq tay0},1); tay0 := mvar numr tay0; % asin's argument l0 := make!-cst!-powerlist TayTemplate tay0; c0 := TayGetCoeff(l0,TayCoeffList tay0); if car u memq '(asec acsc asech acsch) then if null numr c0 then return taysimpasec!*(car u,tay0) else tay := invtaylor tay0 else tay := tay0; tp := TayTemplate tay; l := prune!-coefflist TayCoeffList tay; if null l then return !*tay2q cst!-Taylor!*(simp!* {car u,0},tp); if is!-neg!-pl TayCfPl car l then return taysimpasin!*(car u,tay); tay2 := multtaylor(tay,tay); if car u memq '(asin acos acsc asec) then tay2 := negtaylor tay2; tay2 := addtaylor( cst!-Taylor!*( !*f2q(if car u memq '(acosh asech) then -1 else 1), tp), tay2); if Taylor!*!-zerop tay2 then Taylor!-error!*('branch!-point,car u) else if null numr TayGetCoeff(l0,TayCoeffList tay2) then return taysimpasin!*(car u,tay); tay2 := invtaylor expttayrat(tay2,1 ./ 2); if car u eq '(acos asec) then tay2 := negtaylor tay2; l := for each krnl in TayVars tay collect inttaylorwrttayvar( multtaylor(difftaylorwrttayvar(tay,krnl),tay2), krnl); v := TayCoeffList!-union for each pp in l collect TayCoeffList cdr pp; singlist := nil ./ 1; for each pp in l do if car pp then singlist := addsq(singlist,car pp); % % special treatment for zeroth coefficient % c0 := simp {car u,mk!*sq c0}; v := TayMakeCoeff(l0,c0) . v; tay := make!-Taylor!*( v, tp, if !*taylorkeeporiginal and TayOrig tay then simp {car u,mk!*sq TayOrig tay} else nil, TayFlags tay); if null numr singlist then return !*tay2q tay; if !*taylorkeeporiginal and TayOrig tay then set!-TayOrig(tay,subtrsq(TayOrig tay,singlist)); return addsq(singlist,!*tay2q tay) end; symbolic procedure taysimpasec!*(fn,tay); begin scalar result,tay1,tay2,i1; i1 := simp 'i; if fn memq '(asin acsc) then tay := multtaylorsq(tay,i1); tay1 := cst!-Taylor!*(1 ./ 1,TayTemplate tay); tay2 := multtaylor(tay,tay); if fn memq '(asec asech) then tay2 := negtaylor tay2; result := taysimplog {'log, addtaylor( expttayrat(addtaylor(tay2,tay1),1 ./ 2), tay1)}; tay1 := taysimplog {'log,tay}; if fn memq '(asin acos asec acsc) then <<result := taysimpsq negsq multsq(result,i1); result := addtaylor!-as!-sq(result,multsq(i1,tay1))>> else result := addtaylor!-as!-sq(result, negsq taysimplog {'log,tay}); return taysimpsq!* result end; symbolic procedure taysimpasin!*(fn,tay); begin scalar result,tay1; if fn memq '(asin acsc) then tay := multtaylorsq(tay,simp 'i); tay1 := cst!-Taylor!*( (if fn memq '(asin asinh acsc acsch) then 1 else -1) ./ 1, TayTemplate tay); result := taysimplog {'log, addtaylor( expttayrat(addtaylor(multtaylor(tay,tay), tay1), 1 ./ 2), tay)}; if fn memq '(asin acos asec acsc) then result := quotsq(result,simp 'i); return taysimpsq!* result end; put('asin,'taylorsimpfn,'taysimpasin); put('acos,'taylorsimpfn,'taysimpasin); put('asec,'taylorsimpfn,'taysimpasin); put('acsc,'taylorsimpfn,'taysimpasin); put('asinh,'taylorsimpfn,'taysimpasin); put('acosh,'taylorsimpfn,'taysimpasin); put('asech,'taylorsimpfn,'taysimpasin); put('acsch,'taylorsimpfn,'taysimpasin); symbolic procedure taysimpatan u; if not (car u memq '(atan acot atanh acoth)) or cddr u then confusion 'taysimpatan else begin scalar l,l0,c0,v,tay,tay2,tp,singlist; tay := taysimpsq simp!* cadr u; if not Taylor!-kernel!-sq!-p tay then return mksq({car u,mk!*sq tay},1); tay := mvar numr tay; % atan's argument tp := TayTemplate tay; l0 := make!-cst!-powerlist tp; l := prune!-coefflist TayCoeffList tay; if null l then return !*tay2q cst!-Taylor!*(simp!* {car u,0},tp); if is!-neg!-pl TayCfPl car l then return taysimpatan!*(car u,tay); c0 := get!-cst!-coeff tay; if car u memq '(atan acot) then c0 := subs2!* multsq(c0,simp 'i); if c0 = (1 ./ 1) or c0 = (-1 ./ 1) then return taysimpatan!*(car u,tay); tay2 := multtaylor(tay,tay); if car u memq '(atanh acoth) then tay2 := negtaylor tay2; tay2 := invtaylor addtaylor(cst!-Taylor!*(1 ./ 1,tp),tay2); if car u eq 'acot then tay2 := negtaylor tay2; l := for each krnl in TayVars tay collect inttaylorwrttayvar( multtaylor(difftaylorwrttayvar(tay,krnl),tay2), krnl); v := TayCoeffList!-union for each pp in l collect TayCoeffList cdr pp; singlist := nil ./ 1; for each pp in l do if car pp then singlist := addsq(singlist,car pp); % % special treatment for zeroth coefficient % c0 := simp {car u, mk!*sq TayGetCoeff(l0,TayCoeffList tay)}; v := TayMakeCoeff (l0,c0) . v; tay := make!-Taylor!*( v, tp, if !*taylorkeeporiginal and TayOrig tay then simp {car u,mk!*sq TayOrig tay} else nil, TayFlags tay); if null numr singlist then return !*tay2q tay; if !*taylorkeeporiginal and TayOrig tay then set!-TayOrig(tay,subtrsq(TayOrig tay,singlist)); return addsq(singlist,!*tay2q tay) end; symbolic procedure taysimpatan!*(fn,tay); begin scalar result,tay1; if fn memq '(atan acot) then tay := multtaylorsq(tay,simp 'i); tay1 := cst!-Taylor!*(1 ./ 1,TayTemplate tay); tay := quottaylor(addtaylor(tay1,tay), addtaylor(tay1,negtaylor tay)); result := multsq(taysimplog {'log,tay},1 ./ 2); if fn eq 'atan then result := quotsq(result,simp 'i) else if fn eq 'acot then result := multsq(result,simp 'i); return taysimpsq!* result end; put('atan,'taylorsimpfn,'taysimpatan); put('acot,'taylorsimpfn,'taysimpatan); put('atanh,'taylorsimpfn,'taysimpatan); put('acoth,'taylorsimpfn,'taysimpatan); comment For the logarithm and exponential we use the extension of an algorithm quoted by Knuth who shows how to do this for series in one expansion variable. We extended this to the case of several variables which is straightforward except for one point, see below. Knuth's algorithm works as follows: Assume you want to compute the series W(x) where W(x) = log V(x) Differentiation of this equation gives V'(x) W'(x) = ----- , or V'(x) = W'(x)V(x) . V(x) You make now the ansatz ----- \ n W(x) = > W x , / n ----- substitute this into the above equation and compare powers of x. This yields the recursion formula n-1 V ----- n 1 \ W = ---- - ------ > m W V . n V n V / m n-m 0 0 ----- m=0 The first coefficient must be calculated directly, it is W = log V . 0 0 To use this for series in more than one variable you have to calculate all partial derivatives: n and m refer then to the corresponding component of the multi index. Looking closely one finds that there is an ambiguity: the same coefficient can be calculated using any of the partial derivatives. The only restriction is that the corresponding component of the multi index must not be zero, since we have to divide by it. We resolve this ambiguity by simply taking the first nonzero component. The case of the exponential is nearly the same: differentiation gives W'(x) = V'(x) W(x) , from which we derive the recursion formula n-1 ----- \ n-m W = > --- W V , W = exp V . n / m m n-m 0 0 ----- m=0 ; symbolic procedure taysimplog u; % % Special Taylor expansion function for logarithms % if not (car u eq 'log) or cddr u then confusion 'taysimplog else Taylor!: begin scalar a0,clist,coefflis,il,l0,l,tay,tp,csing,singterm; u := simplogi cadr u; if not kernp u then return taysimpsq u; u := mvar numr u; if not (car u eq 'log) then confusion 'taysimplog; u := taysimpsq simp!* cadr u; if not Taylor!-kernel!-sq!-p u then return mksq({'log,mk!*sq u},1); tay := mvar numr u; tp := TayTemplate tay; l0 := make!-cst!-powerlist tp; % % The following relies on the standard ordering of the % TayCoeffList. % l := prune!-coefflist TayCoeffList tay; if null l then Taylor!-error!*('not!-a!-unit,'taysimplog); % % The assumption at this point is that the first term is the one % with lowest degree, i.e. dividing by this term yields a series % which starts with a constant term. % if TayCfPl car l neq l0 then <<csing := TayCfPl car l; l := for each pp in l collect begin scalar newpl; newpl := subtr!-degrees(TayCfPl pp,csing); if is!-neg!-pl newpl then Taylor!-error('not!-a!-unit,'taysimplog) else return TayMakeCoeff(newpl,TayCfSQ pp); end; tp := addto!-all!-TayTpElOrders( tp, for each nl in csing collect - get!-degree nl); singterm := simp!* retimes preptaycoeff(csing,tp); if !:minusp lnc numr TayCfSq car l then <<singterm := negsq singterm; l := negtaylor1 l>>>>; a0 := TayGetCoeff(l0,l); clist := {TayMakeCoeff(l0,simplogi mk!*sq a0)}; il := if null l then nlist(1,length tp) else smallest!-increment l; coefflis := makecoeffs0(tp,TpNextList tp,il); if not null coefflis then for each cc in cdr coefflis do begin scalar s,pos,pp,n,n1; s := nil ./ 1; pos := find!-non!-zero cc; n := nth(nth(cc,car pos),cdr pos); pp := makecoeffpairs(l0,cc,TayCfPl car l,il); for each p in pp do << n1 := nth(nth(car p,car pos),cdr pos); s := addsq(s, multsq(!*TayExp2q n1, multsq(TayGetCoeff(car p,clist), TayGetCoeff(cdr p,l))))>>; % for each p in pp addsq % multsq(!*TayExp2q nth(nth(car p,car pos),cdr pos), % multsq(TayGetCoeff(car p,clist), % TayGetCoeff(cdr p,l))); s := subtrsq(TayGetCoeff(cc,l),quotsq(s,!*TayExp2q n)); if not null numr s then clist := TayMakeCoeff(cc,quotsq(s,a0)) . clist; end; tay := make!-Taylor!*( reversip clist, tp, if !*taylorkeeporiginal and TayOrig tay then simplogi mk!*sq TayOrig tay else nil, TayFlags tay); if null csing then return !*tay2q tay; singterm := simplogsq singterm; if Taylor!*!-zerop tay then return singterm; if !*taylorkeeporiginal and TayOrig tay then set!-TayOrig(tay,subtrsq(TayOrig tay,singterm)); return addsq(singterm,!*tay2q tay) end; put('log,'taylorsimpfn,'taysimplog); symbolic procedure taysimpexp u; % % Special Taylor expansion function for exponentials % if not (car u eq 'exp) or cddr u then confusion 'taysimpexp else Taylor!: begin scalar a0,clist,coefflis,il,l0,l,lm,lp,tay,tp; u := taysimpsq simp!* cadr u; if not Taylor!-kernel!-sq!-p u then return mksq ({'exp,mk!*sq u},1); tay := mvar numr u; tp := TayTemplate tay; l0 := make!-cst!-powerlist tp; % % The following relies on the standard ordering of the % TayCoeffList. % l := prune!-coefflist TayCoeffList tay; if null l then return !*tay2q cst!-Taylor!*(1 ./ 1,tp); for each pp in l do if is!-neg!-pl TayCfPl pp then lm := pp . lm else if not null numr TayCfSq pp then lp := pp . lp; lm := reversip lm; l := reversip lp; if lm then lm := simp!* {'exp, preptaylor!* make!-Taylor!*(lm,tp,nil,nil)}; if null l then return lm; a0 := TayGetCoeff(l0,l); clist := {TayMakeCoeff(l0,simp!* {'exp,mk!*sq a0})}; il := smallest!-increment l; coefflis := makecoeffs0(tp,TpNextList tp,il); if not null coefflis then for each cc in cdr coefflis do begin scalar s,pos,pp,n,n1; s := nil ./ 1; pos := find!-non!-zero cc; n := nth(nth(cc,car pos),cdr pos); pp := makecoeffpairs(l0,cc,l0,il); for each p in pp do << n1 := nth(nth(car p,car pos),cdr pos); s := addsq(s, multsq(!*TayExp2q(n - n1), multsq(TayGetCoeff(car p,clist), TayGetCoeff(cdr p,l))))>>; s := subs2!* quotsq(s,!*TayExp2q n); if not null numr s then clist := TayMakeCoeff(cc,s) . clist end; clist := reversip clist; u := !*tay2q make!-Taylor!*( clist, tp, if !*taylorkeeporiginal and TayOrig tay then simp {'exp,mk!*sq TayOrig tay} else nil, TayFlags tay); return if null lm then u else multsq(u,lm) end; put('exp,'taylorsimpfn,'taysimpexp); comment The algorithm for the trigonometric functions is also derived from their differential equation. The simplest case is that of tangent whose equation is 2 tan'(x) = 1 + tan (x) . (*) For the others we have 2 cot'(x) = - (1 + cot (x)), 2 tanh'(x) = 1 - tanh (x), 2 coth'(x) = 1 - coth (x) . Let T(x) be a Taylor series, ----- \ N T(x) = > a x / N ----- N=0 Now, let ----- \ N T1(x) = tan T(x) = > b x / N ----- N=0 from which we immediately deduce that b = tan a . 0 0 From (*) we get 2 T1'(x) = (1 + T1(x) ) T'(x) , or, written in terms of the series: Inserting this into (*) we get ----- / / ----- \ 2 \ ----- \ N-1 | | \ N | | \ L > N b x = | 1 + | > b x | | > L a x / N | | / N | | / L ----- \ \ ----- / / ----- N=1 N=0 L=1 We perform the square on the r.h.s. using Cauchy's rule and obtain: ----- \ N-1 > N b x = / N ----- N=1 N / ----- ----- \ ----- | \ \ N | \ L | 1 + > > b b x | > L a x . | / / N-M M | / L \ ----- ----- / ----- N=0 M=0 L=1 Expanding this once again with Cauchy's product rule we get ----- \ N-1 > N b x = / N ----- N=1 L-1 N ----- / ----- ----- \ \ L-1 | \ \ | > x | L a + > > b b (L-N) a | . / | L / / N-M M L-N | ----- \ ----- ----- / L=1 N=0 M=0 From this we immediately deduce the recursion relation L-1 N ----- ----- \ L-N \ b = a + > ----- a > b b . (**) L L / L L-N / N-M M ----- ----- N=0 M=0 This relation is easily generalized to the case of a series in more than one variable, where the same comments apply as in the case of log and exp above. For the hyperbolic tangent the relation is nearly the same. Since its differential equation has a `-' where that of tangent has a `+' we simply have to do the same substitution in the relation (**). For the cotangent we get an additional overall minus sign. There is one additional problem to be handled: if the constant term of T(x), i.e. T(x0) is a pole of tangent. This can be solved quite easily: for a small quantity TAYEPS we calculate Te(x) = T(x) + TAYEPS . and perform the above calculation for Te(x). Then we recover the desired result via the relation tan T(x) = tan (Te(x) - TAYEPS) tan Te(x) - tan TAYEPS = ---------------------------- . 1 + tan Te(x) * tan TAYEPS For the other functions we have similar relations: tanh T(x) = tanh (Te(x) - TAYEPS) tanh Te(x) - tanh TAYEPS = ------------------------------ , 1 - tanh Te(x) * tanh TAYEPS cot T(x) = cot (Te(x) - TAYEPS) 1 + cot Te(x) * cot TAYEPS = ---------------------------- , cot TAYEPS - cot Te(x) coth T(x) = coth (Te(x) - TAYEPS) 1 - coth Te(x) * coth TAYEPS = ------------------------------ . coth Te(x) - coth TAYEPS We know that this result is independent of TAYEPS, and we can thus evaluate it for any value of TAYEPS. Since we further know that T(x0) is a pole of the function in question, we can express tan (T(x0) + TAYEPS) as 1 - ------------ , tan TAYEPS and similarly all the other possible expressions involving TAYEPS can be written in terms of tan TAYEPS and tanh TAYEPS, respectively. This makes it possible to just substitute any occurrence of tan TAYEPS or tanh TAYEPS by zero, which greatly simplifies the final calculation ; !*!*taylor!-epsilon!*!* := compress '(t a y e p s); symbolic procedure taysimptan u; % % Special Taylor expansion function for circular and hyperbolic % tangent and cotangent % if not (car u memq '(tan cot tanh coth)) or cddr u then confusion 'taysimptan else Taylor!: begin scalar a,a0,clist,coefflis,il,l0,l,poleflag,tay,tp; tay := taysimpsq simp!* cadr u; if not Taylor!-kernel!-sq!-p tay then return mksq({car u,mk!*sq tay},1); tay := mvar numr tay; tp := TayTemplate tay; l0 := make!-cst!-powerlist tp; % % First we get rid of possible zero coefficients. % l := prune!-coefflist TayCoeffList tay; % if null l then return !*tay2q cst!-Taylor!*(simp!* {car u,0},tp); % % The following relies on the standard ordering of the % TayCoeffList. % if not null l and is!-neg!-pl TayCfPl car l then Taylor!-error('essential!-singularity,car u); a0 := TayGetCoeff(l0,l); il := if null l then nlist(1,length tp) else smallest!-increment l; % %%% handle poles of function % a := quotsq(a0,simp 'pi); if car u memq '(tanh coth) then a := subs2!* multsq(a,simp 'i); if car u memq '(tan tanh) and denr(a := multsq(a,simp '2)) = 1 and fixp (a := numr a) and not evenp a or car u memq '(cot coth) and denr a = 1 and (null (a := numr a) or fixp a) then << % % OK, now we are at a pole, so we add a small quantity, compute % the series and use the addition formulas to get the real % result. % poleflag := t; a0 := if car u eq 'tan then negsq invsq simp!* {'tan,!*!*taylor!-epsilon!*!*} else if car u eq 'tanh then invsq simp!* {'tanh,!*!*taylor!-epsilon!*!*} else if car u eq 'cot then invsq simp!* {'tan,!*!*taylor!-epsilon!*!*} else invsq simp!* {'tanh,!*!*taylor!-epsilon!*!*}; clist := {TayMakeCoeff(l0,a0)}; >> else clist := {TayMakeCoeff(l0,simp!* {car u,mk!*sq a0})}; % coefflis := makecoeffs0(tp,TpNextList tp,il); if not null coefflis then for each cc in cdr coefflis do begin scalar cf,s,pos,x,y,n,n1; s := nil ./ 1; pos := find!-non!-zero cc; n := nth(nth(cc,car pos),cdr pos); for each p in makecoeffpairs(l0,cc,l0,il) do << x := reversip makecoeffpairs1(l0,car p,l0,il); y := nil ./ 1; for each z in x do y := addsq(y,multsq(TayGetCoeff(car z,clist), TayGetCoeff(cdr z,clist))); n1 := nth(nth(car p,car pos),cdr pos); s := addsq(s, multsq(!*TayExp2q(n - n1), multsq(y,TayGetCoeff(cdr p,l))))>>; cf := quotsq(s,!*TayExp2q n); if car u memq '(tanh coth) then cf := negsq cf; cf := addsq(TayGetCoeff(cc,l),cf); if null numr cf then return; % short cut for efficiency if car u eq 'cot then cf := negsq cf; clist := TayMakeCoeff(cc,cf) . clist end; a := make!-Taylor!*(reversip clist,tp,nil,nil); % % Construct ``real'' series in case of pole % if poleflag then begin scalar x; x := if car u eq 'cot then cst!-Taylor!*( invsq simp {'tan,!*!*taylor!-epsilon!*!*},tp) else if car u eq 'coth then cst!-Taylor!*( invsq simp {'tanh,!*!*taylor!-epsilon!*!*},tp) else cst!-Taylor!*( simp {car u,!*!*taylor!-epsilon!*!*},tp); if car u eq 'tan then a := quottaylor(addtaylor(a,negtaylor x), addtaylor(cst!-taylor!*(1 ./ 1,tp), multtaylor(a,x))) else if car u eq 'cot then a := quottaylor(addtaylor(multtaylor(a,x), cst!-taylor!*(1 ./ 1,tp)), addtaylor(x,negtaylor a)) else if car u eq 'tanh then a := quottaylor(addtaylor(a,negtaylor x), addtaylor(cst!-taylor!*(1 ./ 1,tp), negtaylor multtaylor(a,x))) else if car u eq 'coth then a := quottaylor(addtaylor(multtaylor(a,x), cst!-taylor!*(-1 ./ 1,tp)), addtaylor(x,negtaylor a)); if not (a freeof !*!*taylor!-epsilon!*!*) then set!-TayCoeffList(a, for each pp in TayCoeffList a collect TayMakeCoeff(TayCfPl pp, subsq(TayCfSq pp, {!*!*taylor!-epsilon!*!* . 0}))); end; % if !*taylorkeeporiginal and TayOrig tay then set!-TayOrig(a,simp {car u,mk!*sq TayOrig tay}); set!-Tayflags(a,TayFlags tay); return !*tay2q a end; put('tan,'taylorsimpfn,'taysimptan); put('cot,'taylorsimpfn,'taysimptan); put('tanh,'taylorsimpfn,'taysimptan); put('coth,'taylorsimpfn,'taysimptan); comment For the circular sine and cosine and their reciprocals we calculate the exponential and use it via the formulas exp(+I*x) - exp(-I*x) sin x = ----------------------- , 2 exp(+I*x) + exp(-I*x) cos x = ----------------------- , 2 etc. To avoid clumsy expressions we separate the constant term of the argument, T(x) = a0 + T1(x), and use the addition theorems which give 1 sin T(x) = - exp(+I*T1(x)) * (sin a0 - I*cos a0) + 2 1 - exp(-I*T1(x)) * (sin a0 + I*cos a0) , 2 1 cos T(x) = - exp(+I*T1(x)) * (cos a0 + I*sin a0) + 2 1 - exp(-I*T1(x)) * (cos a0 - I*sin a0) . 2 ; symbolic procedure taysimpsin u; % % Special Taylor expansion function for circular sine, cosine, % and their reciprocals % if not (car u memq '(sin cos sec csc)) or cddr u then confusion 'taysimpsin else Taylor!: begin scalar l,tay,result,tp,i1,l0,a0,a1,a2; tay := taysimpsq simp!* cadr u; if not Taylor!-kernel!-sq!-p tay then return mksq({car u,mk!*sq tay},1); tay := mvar numr tay; tp := TayTemplate tay; l0 := make!-cst!-powerlist tp; l := prune!-coefflist TayCoeffList tay; % if null l then return !*tay2q cst!-Taylor!*(simp!* {car u,0},tp); % if is!-neg!-pl TayCfPl car l % then Taylor!-error('essential!-singularity,car u); a0 := TayGetCoeff(l0,l); % % make constant term to 0 % i1 := simp 'i; if not null numr a0 then tay := addtaylor(tay,cst!-Taylor!*(negsq a0,tp)); result := taysimpexp{'exp,multtaylor(tay,cst!-Taylor!*(i1,tp))}; a1 := simp!* {'sin,mk!*sq a0} . simp!* {'cos,mk!*sq a0}; if car u memq '(sin csc) then << a2 := addsq(car a1,multsq(i1,cdr a1)); a1 := addsq(car a1,negsq multsq(i1,cdr a1)); >> else << a2 := addsq(cdr a1,negsq multsq(i1,car a1)); a1 := addsq(cdr a1,multsq(i1,car a1)); >>; result := multsq(addsq(multsq(result,a1), multsq(taysimpsq!* invsq result,a2)), 1 ./ 2); if car u memq '(sec csc) then result := invsq result; return taysimpsq!* result end; put('sin,'taylorsimpfn,'taysimpsin); put('cos,'taylorsimpfn,'taysimpsin); put('sec,'taylorsimpfn,'taysimpsin); put('csc,'taylorsimpfn,'taysimpsin); comment For the hyperbolic sine and cosine and their reciprocals we calculate the exponential and use it via the formulas exp(+x) - exp(-x) sinh x = ------------------- , 2 exp(+x) + exp(-x) cosh x = ------------------- , 2 etc. To avoid clumsy expressions we separate the constant term of the argument, T(x) = a0 + T1(x), and use the addition theorems which give 1 sinh T(x) = - exp(+T1(x)) * (sinh a0 + cosh a0) + 2 1 - exp(-T1(x)) * (sinh a0 - cosh a0) , 2 1 cosh T(x) = - exp(+T1(x)) * (cosh a0 + sinh a0) + 2 1 - exp(-T1(x)) * (cosh a0 - sinh a0) . 2 ; symbolic procedure taysimpsinh u; % % Special Taylor expansion function for circular sine, cosine, % and their reciprocals % if not (car u memq '(sinh cosh sech csch)) or cddr u then confusion 'taysimpsin else Taylor!: begin scalar l,tay,result,tp,l0,a0,a1,a2; tay := taysimpsq simp!* cadr u; if not Taylor!-kernel!-sq!-p tay then return mksq({car u,mk!*sq tay},1); tay := mvar numr tay; tp := TayTemplate tay; l0 := make!-cst!-powerlist tp; l := prune!-coefflist TayCoeffList tay; % if null l then return !*tay2q cst!-Taylor!*(simp!* {car u,0},tp); % if is!-neg!-pl TayCfPl car l % then Taylor!-error('essential!-singularity,car u); a0 := TayGetCoeff(l0,l); % % make constant term to 0 % if not null numr a0 then tay := addtaylor(tay,cst!-Taylor!*(negsq a0,tp)); result := taysimpexp{'exp,tay}; a1 := simp!* {'sinh,mk!*sq a0} . simp!* {'cosh,mk!*sq a0}; if car u memq '(sinh csch) then << a2 := addsq(car a1,cdr a1); a1 := subtrsq(car a1,cdr a1); >> else << a2 := addsq(cdr a1,car a1); a1 := subtrsq(cdr a1,car a1); >>; result := multsq(addsq(multsq(result,a2), multsq(taysimpsq!* invsq result,a1)), 1 ./ 2); if car u memq '(sech csch) then result := invsq result; return taysimpsq!* result end; put('sinh, 'taylorsimpfn, 'taysimpsinh); put('cosh, 'taylorsimpfn, 'taysimpsinh); put('sech, 'taylorsimpfn, 'taysimpsinh); put('csch, 'taylorsimpfn, 'taysimpsinh); comment Support for the integration of Taylor kernels. Unfortunately, with the current interface, only Taylor kernels on toplevel can be treated successfully. The way it is down means stretching certain interfaces beyond what they were designed for...but it works! First we add a rule that replaces a call to INT with a Taylor kernel as argument by a call to TAYLORINT, then we define REVALTAYLORINT as simplification function for that; algebraic let { int(symbolic algebraic(taylor!*(~x,~y,~z,~w)),~u) => taylorint(~x,~y,~z,~w,~u), int(log(~u)^~~n * symbolic algebraic(taylor!*(~x,~y,~z,~w)),~u) => log(u)^n * int(symbolic('(taylor!* x y z w)),u) - int(log(u)^(n-1) * taylorcombine(int(symbolic('(taylor!* x y z w)),u)/u),u), int(~x,~y) => taylorint1(~x,~y) when not symbolic algebraic(~x freeof 'Taylor!*) }; put('taylorint1, 'psopfn, 'revaltaylorint1); symbolic procedure revaltaylorint1 x; begin scalar u,v; u := car x; v := cadr x; if Taylor!*p u then return revaltaylorint append(cdr u,{v}); u := reval taylorcombine u; if Taylor!*p u then return revaltaylorint append(cdr u,{v}); if not atom u then if car u memq '(plus minus difference) then return reval (car u . for each term in cdr u collect {'int,term,v}); lprim "Converting Taylor kernels to standard representation"; return aeval {'int,taylortostandard car x,v} end; put('taylorint, 'psopfn, 'revaltaylorint); symbolic procedure revaltaylorint u; begin scalar taycfl, taytp, tayorg, tayflgs, var; taycfl := car u; taytp := cadr u; tayorg := caddr u; tayflgs := cadddr u; var := car cddddr u; if null (var member TayTpVars taytp) then return mk!*sq !*tay2q make!-Taylor!*( for each pp in taycfl collect TayCfPl pp . simp!* {'int,mk!*sq TayCfSq pp,var}, taytp, if not null tayorg then simp!* {'int,mk!*sq tayorg,var} else nil, nil) else return mk!*sq ((if null car result then !*tay2q cdr result else addsq(car result,!*tay2q cdr result)) where result := inttaylorwrttayvar1(taycfl,taytp,var)) end; endmodule; end;