Artifact 0258d77603422bbe2ebe6fb967049a07b4e9ddce6f075d13a46be5a1e823fd4b:
- Executable file
r37/packages/int/halfangl.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: 10996) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/int/halfangl.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: 10996) [annotate] [blame] [check-ins using]
module halfangl; % Routines for conversion to half angle tangents. % Author: Steve Harrington. % Modifications by: John P. Fitch. fluid '(!*gcd); exports halfangle,untan; symbolic procedure transform(u,x); % Transform the SQ U to remove the 'bad' functions sin, cos, cot etc. % in favor of half angles. % Do this with regard to cases like sin(x)*tan(x/2), so attempt to % limit times we use half angles. begin scalar zl,tnarg,substlist; zl := zlist; while car(tnarg := tan!-function!-in zl) and halfangle!-confusion(zlist,cadar tnarg) do <<zl := cdr tnarg; tnarg := car tnarg; if eqcar(tnarg,'tan) then substlist := (gensym() . tnarg) . substlist else substlist := (gensym() . list('quotient,1,('tan . cdr tnarg))) . substlist; u := subst(caar substlist,tnarg,u)>>; return if substlist then simp sublis(substlist,prepsq halfangle(u,x)) % simp prepsq was added so that 1/(e**x*cos(1/e**x)**2) % for example returns a closed-form result. else simp prepsq halfangle(u,x) end; symbolic procedure tan!-function!-in zz; % Look at zlist for tangents or cotangents. <<while zz and not eqcar(car zz,'tan) and not eqcar(car zz,'cot) do zz := cdr zz; if null zz then nil . nil else zz>>; symbolic procedure halfangle!-confusion(zz,tnarg); % Is there a function in the zlist with twice the tangent argument? <<while zz and (atom car zz or not(tnarg = prepsq simp list('quotient,cadar zz,2))) do zz := cdr zz; zz>>; symbolic procedure quotqq(u1,v1); multsq(u1,invsq(v1)); symbolic procedure !*subtrq(u1,v1); addsq(u1, negsq(v1)); symbolic procedure !*int2qm(u1); if u1=0 then nil . 1 else u1 . 1; symbolic procedure halfangle(r,x); % Top level procedure for converting; % R is a rational expression to be converted, % X the integration variable. % A rational expression is returned. quotqq(hfaglf(numr(r),x), hfaglf(denr(r),x)); symbolic procedure hfaglf(p,x); % Converting polynomials, a rational expression is returned. if domainp(p) then !*f2q(p) else subs2q addsq(multsq(exptsq(hfaglk(mvar(p),x), ldeg(p)), hfaglf(lc(p),x)), hfaglf(red(p),x)); symbolic procedure hfaglk(k,x); % Converting kernels, a rational expression is returned. begin scalar kt; if atom k or not member(x,flatten(cdr(k))) then return !*k2q k; k := car(k) . hfaglargs(cdr(k), x); if cadr k eq 'pi then return !*k2q k; % Don't consider tan(pi/2). kt := simp list('tan, list('quotient, cadr(k), 2)); return if car(k) = 'sin then quotqq(multsq(!*int2qm(2),kt), addsq(!*int2qm(1), exptsq(kt,2))) else if car(k) = 'cos then quotqq(!*subtrq(!*int2qm(1),exptsq(kt,2)),addsq(!*int2qm(1), exptsq(kt,2))) else if car(k) = 'tan then quotqq(multsq(!*int2qm(2),kt), !*subtrq(!*int2qm(1), exptsq(kt,2))) else if car(k) = 'cot then quotqq(!*subtrq(!*int2qm(1), exptsq(kt,2)),multsq(!*int2qm(2),kt)) else if car(k) = 'sec then quotqq(addsq(!*int2qm(1), exptsq(kt,2)), !*subtrq(!*int2qm(1),exptsq(kt,2))) else if car(k) = 'csc then quotqq(addsq(!*int2qm(1),exptsq(kt,2)), %%% !*subtrq(!*int2qm(1),exptsq(kt,2))) % FJW - was identical to sec!!! multsq(!*int2qm(2),kt)) else if car(k) = 'sinh then quotqq(!*subtrq(!*p2q mksp('expt.('e. cdr k),2), !*int2qm(1)), multsq(!*int2qm(2), !*p2q mksp('expt . ('e . cdr(k)),1))) else if car(k) = 'cosh then quotqq(addsq(!*p2q mksp('expt.('e. cdr k),2), !*int2qm(1)), multsq(!*int2qm(2), !*p2q mksp('expt . ('e . cdr(k)),1))) else if car(k) = 'tanh then quotqq(!*subtrq(!*p2q mksp('expt.('e. cdr k),2), !*int2qm(1)), addsq(!*p2q mksp ('expt.('e.cdr(k)),2), !*int2qm(1))) else if car(k) = 'coth then quotqq(addsq(!*p2q mksp('expt.('e.cdr(k)),2), !*int2qm(1)), !*subtrq(!*p2q mksp('expt.('e . cdr k),2),!*int2qm(1))) else if car(k) = 'acot then !*p2q mksp(list('atan, list('quotient, 1, cadr k)),1) else !*k2q(k); % additional transformation might be added here. end; symbolic procedure hfaglargs(l,x); % Conversion of argument list. if null l then nil else prepsq(hfaglk(car(l),x)) . hfaglargs(cdr(l),x); symbolic procedure untanf x; % This should be done by a table. % We turn off gcd to avoid unnecessary gcd calculations, as suggested % by Rainer Schoepf. begin scalar !*gcd,y,z,w; if domainp x then return x . 1; y := mvar x; if eqcar(y,'int) then error1(); % assume all is hopeless. z := ldeg x; w := 1 . 1; y := if atom y then !*k2q y else if car y eq 'tan then begin scalar yy; %% printc "Recursive tan"; printc cadr y; yy := prepsq untan simp cadr y . nil; %% princ "==> "; printc yy; if evenp z then <<z := z/2; return simp list('quotient, list('plus, list('minus, list('cos, 'times . (2 . yy))), 1),list('plus, list('cos, 'times . (2 . yy)), 1))>> else if z=1 then return simp list('quotient, list('plus, list('minus, list('cos, 'times . (2 . yy))), 1),list('sin, 'times . (2 . yy))) else <<z := (z - 1)/2; w := simp list('quotient, list('plus, list('minus, list('cos, 'times . (2 . yy))), 1),list('sin, 'times . (2 . yy))); return simp list('quotient, list('plus, list('minus, list('cos, 'times . (2 . yy))), 1),list('plus, list('cos, 'times . (2 . yy)), 1)) >> end else simp y; return addsq(multsq(multsq(exptsq(y,z),untanf lc x),w), untanf red x) end; % symbolic procedure untanlist(y); % if null y then nil % else (prepsq (untan(simp car y)) . untanlist(cdr y)); symbolic procedure untan(x); % Expects x to be canonical quotient. begin scalar y; y:=cossqchk sinsqrdchk multsq(untanf(numr x), invsq untanf(denr x)); return if length flatten y>length flatten x then x else y end; symbolic procedure sinsqrdchk(x); multsq(sinsqchkf(numr x), invsq sinsqchkf(denr x)); symbolic procedure sinsqchkf(x); begin scalar y,z,w; if domainp x then return x . 1; y := mvar x; z := ldeg x; w := 1 . 1; y := if eqcar(y,'sin) then if evenp z then <<z := quotient(z,2); simp list('plus,1,list('minus, list('expt,('cos . cdr(y)),2)))>> else if z = 1 then !*k2q y else << z := quotient(difference(z,1),2); w := !*k2q y; simp list('plus,1,list('minus, list('expt,('cos . cdr(y)),2)))>> else !*k2q y; return addsq(multsq(multsq(exptsq(y,z),sinsqchkf(lc x)),w), sinsqchkf(red x)); end; symbolic procedure cossqchkf(x); begin scalar y,z,w,x1,x2; if domainp x then return x . 1; y := mvar x; z := ldeg x; w := 1 . 1; x1 := cossqchkf(lc x); x2 := cossqchkf(red x); x := addsq(multsq(!*p2q lpow x,x1),x2); y := if eqcar(y,'cos) then if evenp z then <<z := quotient(z,2); simp list('plus,1,list('minus, list('expt,('sin . cdr(y)),2)))>> else if z = 1 then !*k2q y else << z := quotient(difference(z,1),2); w := !*k2q y; simp list('plus,1,list('minus, list('expt,('sin . cdr(y)),2)))>> else !*k2q y; y := addsq(multsq(multsq(exptsq(y,z),w),x1),x2); return if length(y) > length(x) then x else y; end; symbolic procedure cossqchk(x); begin scalar !*gcd; !*gcd := t; return multsq(cossqchkf(numr x), invsq cossqchkf(denr x)) end; symbolic procedure lrootchk(l,x); % Checks each member of list l for a root. if null l then nil else krootchk(car l, x) or lrootchk(cdr l, x); symbolic procedure krootchk(f,x); % Checks a kernel to see if it is a root. if atom f then nil else if car(f) = 'sqrt and member(x, flatten cdr f) then t else if car(f) = 'expt and not atom caddr(f) and caaddr(f) = 'quotient and member(x, flatten cadr f) then t else lrootchk(cdr f, x); symbolic procedure rootchk1p(f,x); % Checks polynomial for a root. if domainp f then nil else krootchk(mvar f,x) or rootchk1p(lc f,x) or rootchk1p(red f,x); symbolic procedure rootcheckp(f,x); % Checks rational (standard quotient) for a root. rootchk1p(numr f,x) or rootchk1p(denr f,x); endmodule; end;