Artifact e080aac1e3a2339c9084a27514d05c764abcd1cc8bb934e67015c964fa157c9a:
- Executable file
r36/src/trigsimp.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: 17682) [annotate] [blame] [check-ins using] [more...]
module trigsimp; % User controlled simplification % of trigonometric expressions. % Authors: Wolfram Koepf, Andreas Bernig, Herbert Melenk % Version 1.0 April 1995 !#if (memq 'psl lispsystem!*) flag('(trigsmp1),'lap); !#endif create!-package('(trigsimp trigsmp1 trigsmp2),'(contrib misc)); load_package compact; endmodule; module trigsmp1; % Collection of rule sets. algebraic; clearrules(trig_imag_rules); trig_normalize!*:= { cos(~a)^2 => 1 - sin(a)^2 when trig_preference=sin, sin(~a)^2 => 1 - cos(a)^2 when trig_preference=cos, cosh(~a)^2 => 1+sinh(a)^2 when hyp_preference=sinh, sinh(~a)^2 => cosh(a)^2-1 when hyp_preference=cosh }; trig_expand!*:= { % additions theorems sin((~a + ~b)/~~m) => sin(a/m)*cos(b/m) + cos(a/m)*sin(b/m), cos((~a + ~b)/~~m) => cos(a/m)*cos(b/m) - sin(a/m)*sin(b/m), tan((~a+~b)/~~m) => (tan(a/m)+tan(b/m))/(1-tan(a/m)*tan(b/m)), cot((~a+~b)/~~m) =>(cot(a/m)*cot(b/m)-1)/(cot(a/m)+cot(b/m)), sec((~a+~b)/~~m) =>1/(1/(sec(a/m)*sec(b/m))-1/(csc(a/m)*csc(b/m))), csc((~a+~b)/~~m) =>1/(1/(sec(b/m)*csc(a/m))+1/(sec(a/m)*csc(b/m))), tanh((~a+~b)/~~m) => (tanh(a/m)+tanh(b/m))/(1+tanh(a/m)*tanh(b/m)), coth((~a+~b)/~~m) =>(coth(a/m)*coth(b/m)+1)/(coth(a/m)+coth(b/m)), sinh((~a+~b)/~~m) => sinh(a/m)*cosh(b/m) + cosh(a/m)*sinh(b/m), cosh((~a+~b)/~~m) => cosh(a/m)*cosh(b/m) + sinh(a/m)*sinh(b/m), sech((~a+~b)/~~m) =>1/(1/(sech(a/m)*sech(b/m))+1/(csch(a/m)*csch(b/m))), csch((~a+~b)/~~m) =>1/(1/(sech(a/m)*csch(b/m))+1/(sech(b/m)*csch(a/m))), % multiplication theorems sin(~n*~a/~~m) => sin(a/m)*cos((n-1)*a/m) + cos(a/m)*sin((n-1)*a/m) when fixp n and n>1 and n<=15, sin(~n*~a/~~m) =>2*sin(n/2*a/m)*cos(n/2*a/m) when fixp n and mod(n,2)=0 and n>15, sin(~n*~a/~~m) =>sin((n-1)/2*a/m)*cos((n+1)/2*a/m)+ sin((n+1)/2*a/m)*cos((n-1)/2*a/m) when fixp n and mod(n,2)=1 and n>15, cos(~n*~a/~~m) => cos(a/m)*cos((n-1)*a/m) - sin(a/m)*sin((n-1)*a/m) when fixp n and n>1 and n<=15, cos(~n*~a/~~m) => 2*cos(n/2*a/m)**2-1 when fixp n and mod(n,2)=0 and n>15, cos(~n*~a/~~m) => cos((n-1)/2*a/m)*cos((n+1)/2*a/m)- sin((n-1)/2*a/m)*sin((n+1)/2*a/m) when fixp n and mod(n,2)=1 and n>15, sinh(~n*~a/~~m) => sinh(a/m)*cosh((n-1)*a/m)+cosh(a/m)*sinh((n-1)*a/m) when fixp n and n<=15 and n>1, sinh(~n*~a/~~m) => 2*sinh(n/2*a/m)*cosh(n/2*a/m) when fixp n and mod(n,2)=0 and n>15, sinh(~n*~a/~~m) => sinh((n-1)/2*a/m)*cosh((n+1)/2*a/m)+ sinh((n+1)/2*a/m)*cosh((n-1)/2*a/m) when fixp n and mod(n,2)=1 and n>15, cosh(~n*~a/~~m) => cosh(a/m)*cosh((n-1)*a/m) + sinh(a/m)*sinh((n-1)*a/m) when fixp n and n>1 and n<=15, cosh(~n*~a/~~m) => 2*cosh(n/2*a/m)**2-1 when fixp n and mod(n,2)=0 and n>15, cosh(~n*~a/~~m) => cosh((n-1)/2*a/m)*cosh((n+1)/2*a/m)+ sinh((n-1)/2*a/m)*sinh((n+1)/2*a/m) when fixp n and mod(n,2)=1 and n>15, tan(~n*~a/~~m) => (tan(a/m)+tan((n-1)*a/m))/(1-tan(a/m)*tan((n-1)*a/m)) when fixp n and n>1 and n<=15, tan(~n*~a/~~m) => 2*tan(n/2*a/m)/(1-tan(n/2*a/m)**2) when fixp n and mod(n,2)=0 and n>15, tan(~n*~a/~~m) => ( tan((n-1)/2*a/m)+tan((n+1)/2*a/m) )/ (1-tan((n-1)/2*a/m)*tan((n+1)/2*a/m)) when fixp n and mod(n,2)=1 and n>15, tanh(~n*~a/~~m) => (tanh(a/m)+tanh((n-1)*a/m))/(1+tanh(a/m)*tanh((n-1)*a/m)) when fixp n and n>1 and n<=15, tanh(~n*~a/~~m) => 2*tanh(n/2*a/m)/(1+tanh(n/2*a/m)**2) when fixp n and mod(n,2)=0 and n>15, tanh(~n*~a/~~m) => ( tanh((n-1)/2*a/m)+tanh((n+1)/2*a/m) )/ (1+tanh((n-1)/2*a/m)*tanh((n+1)/2*a/m)) when fixp n and mod(n,2)=1 and n>15, cot(~n*~a/~~m) => (cot(a/m)*cot((n-1)*a/m)-1)/(cot(a/m)+cot((n-1)*a/m)) when fixp n and n>1 and n<=15, cot(~n*~a/~~m) => (cot(n/2*a/m)**2-1)/(2cot(n/2*a/m)) when fixp n and mod(n,2)=0 and n>15, cot(~n*~a/~~m) => ( cot((n-1)/2*a/m)*cot((n+1)/2*a/m)-1 ) / (cot((n-1)/2*a/m)+cot((n+1)/2*a/m)) when fixp n and mod(n,2)=1 and n>15, coth(~n*~a/~~m) => (coth(a/m)*coth((n-1)*a/m)+1)/(coth(a/m)+coth((n-1)*a/m)) when fixp n and n>1 and n<=15, coth(~n*~a/~~m) => (coth(n/2*a/m)**2+1)/(2coth(n/2*a/m)) when fixp n and mod(n,2)=0 and n>15, coth(~n*~a/~~m) => ( coth((n-1)/2*a/m)*coth((n+1)/2*a/m)+1 ) / (coth((n-1)/2*a/m)+coth((n+1)/2*a/m)) when fixp n and mod(n,2)=1 and n>15, sec(~n*~a/~~m) => 1/(1/(sec(a/m)*sec((n-1)*a/m))-1/(csc(a/m)*csc((n-1)*a/m))) when fixp n and n>1 and n<=15, sec(~n*~a/~~m) =>1/(1/sec(n/2*a/m)**2-1/csc(n/2*a/m)**2) when fixp n and mod(n,2)=0 and n>15, sec(~n*~a/~~m) =>1/(1/(sec((n-1)/2*a/m)*sec((n+1)/2*a/m))- 1/(csc((n-1)/2*a/m)*csc((n+1)/2*a/m))) when fixp n and mod(n,2)=1 and n>15, csc(~n*~a/~~m) => 1/(1/(sec(a/m)*csc((n-1)*a/m))+1/(csc(a/m)*sec((n-1)*a/m))) when fixp n and n>1 and n<=15, csc(~n*~a/~~m) =>sec(n/2*a/m)*csc(n/2*a/m)/2 when fixp n and mod(n,2)=0, csc(~n*~a/~~m) =>1/(1/(sec((n-1)/2*a/m)*csc((n+1)/2*a/m))+ 1/(csc((n-1)/2*a/m)*sec((n+1)/2*a/m))) when fixp n and mod(n,2)=1 and n>15, sech(~n*~a/~~m) => 1/(1/(sech(a/m)*sech((n-1)*a/m))+1/(csch(a/m)*csch((n-1)*a/m))) when fixp n and n>1 and n<=15, sech(~n*~a/~~m) =>1/(1/sech(n/2*a/m)**2+1/csch(n/2*a/m)**2) when fixp n and mod(n,2)=0 and n>15, sech(~n*~a/~~m) =>1/(1/(sech((n-1)/2*a/m)*sech((n+1)/2*a/m))+ 1/(csch((n-1)/2*a/m)*csch((n+1)/2*a/m))) when fixp n and mod(n,2)=1 and n>15, csch(~n*~a/~~m) => 1/(1/(sech(a/m)*csch((n-1)*a/m))+1/(csch(a/m)*sech((n-1)*a/m))) when fixp n and n>1 and n<=15, csch(~n*~a/~~m) =>sech(n/2*a/m)*csch(n/2*a/m)/2 when fixp n and mod(n,2)=0 and n>15, csch(~n*~a/~~m) =>1/(1/(sech((n-1)/2*a/m)*csch((n+1)/2*a/m))+ 1/(csch((n-1)/2*a/m)*sech((n+1)/2*a/m))) when fixp n and mod(n,2)=1 and n>15 }; trig_combine!*:= { sin(~a)*sin(~b) => 1/2*(cos(a-b) - cos(a+b)), cos(~a)*cos(~b) => 1/2*(cos(a-b) + cos(a+b)), sin(~a)*cos(~b) => 1/2*(sin(a-b) + sin(a+b)), sin(~a)^2 => 1/2*(1-cos(2*a)), cos(~a)^2 => 1/2*(1+cos(2*a)), sinh(~a)*sinh(~b) => 1/2*(cosh(a+b) - cosh(a-b)), cosh(~a)*cosh(~b) => 1/2*(cosh(a-b) + cosh(a+b)), sinh(~a)*cosh(~b) => 1/2*(sinh(a-b) + sinh(a+b)), sinh(~a)^2 => 1/2*(cosh(2*a)-1), cosh(~a)^2 => 1/2*(1+cosh(2*a)) }; trig_standardize!*:= { tan(~a) => sin(a)/cos(a), cot(~a) => cos(a)/sin(a), tanh(~a) => sinh(a)/cosh(a), coth(~a) => cosh(a)/sinh(a), sec(~a) => 1/cos(a), csc(~a) => 1/sin(a), sech(~a) => 1/cosh(a), csch(~a) => 1/sinh(a) } ; trig2exp!*:= { cos(~a) => (e^(i*a) + e^(-i*a))/2, sin(~a) => -i*(e^(i*a) - e^(-i*a))/2, cosh(~a) => (e^(a) + e^(-a))/2, sinh(~a) => (e^(a) - e^(-a))/2 }; exp2trig1!*:= { e**(~x)=>cos(x/i)+i*sin(x/i) }; exp2trig2!*:= { e**(~x)=>1/(cos(x/i)-i*sin(x/i)) }; trig2hyp!*:= { sin(~a)=> -i * sinh(i*a), cos(~a)=> cosh(i*a), tan(~a)=> -i*tanh(i*a), cot(~a)=> i*coth(i*a), sec(~a)=> sech(i*a), csc(~a)=> i*csch(i*a), asin(~a)=> asinh(i*a)/i, acos(~a)=> acosh(a)/i }; hyp2trig!*:= { sinh(~a)=> -i*sin(i*~a), cosh(~a)=> cos(i*~a), acosh(~a)=> i*acos(a), asinh(~a)=> asin(-i*a)*i }; subtan!*:= { sin(~x)=>cos(x)*tan(x) when trig_preference=cos, cos(~x)=>sin(x)/tan(x) when trig_preference=sin, sinh(~x)=>cosh(x)*tanh(x) when hyp_preference=cosh, cosh(x)=>sinh(x)/tanh(x) when hyp_preference=sinh }; endmodule; module trigsmp2; % Executable code. algebraic; procedure behandle(ex); begin scalar p,q,p2,q2; p := num ex; q := den ex; let exp2trig1!*; p2 := p; clearrules exp2trig1!*; let exp2trig2!*; q2 := q; clearrules exp2trig2!*; return p2/q2; end; procedure trigsimp1(f,l); begin scalar u,p,trigpreferencelist,hyppreferencelist, directionlist,modelist,keepalltriglist,err,onlytan; err:=0; if freeof(f,sin) and freeof(f,cos) and freeof(f,cosh) and freeof(f,sinh) and freeof(f,csc) and freeof(f,sec) and freeof(f,csch) and freeof(f,sech) then onlytan:=1 else onlytan:=0; trigpreferencelist:={}; hyppreferencelist:={}; directionlist:={}; modelist:={}; keepalltriglist:={}; while length(l) neq 0 do begin u:=first(l); l:=rest(l); if u=sin or u=cos then trigpreferencelist:=u.trigpreferencelist else if (u=sinh) or (u=cosh) then hyppreferencelist:= u.hyppreferencelist else if (u=expand) or (u=combine) or (u=compact) then directionlist:=u.directionlist else if u=hyp or u=trig or u=expon then modelist:=u.modelist else if (u=keepalltrig) then keepalltriglist:=u.keepalltriglist else <<write u," not possible. Use sin,cos,cosh,sinh, expand,combine,compact,hyp,trig,expon,keepalltrig!";err:=1;>>; end; %Defaults if trigpreferencelist={} then trigpreferencelist:={sin}; if hyppreferencelist={} then hyppreferencelist:={sinh}; if directionlist={} then directionlist:={expand}; %contradictions? if length(trigpreferencelist)>1 then <<write trigpreferencelist, " not possible. Use either sin or cos.";err:=1;>>; if length(hyppreferencelist)>1 then <<write hyppreferencelist, " not possible. Use either sinh or cosh.";err:=1;>>; if length(directionlist)>1 then <<write directionlist, " not possible. Use either expand or combine or compact.";err:=1;>>; if length(modelist)>1 then <<write modelist," not possible. Use either hyp or expon.";err:=1;>>; if err=0 then begin %application if first(trigpreferencelist)=sin then trig_preference:=sin; if first(trigpreferencelist)=cos then trig_preference:=cos; if first(hyppreferencelist)=sinh then hyp_preference:=sinh; if first(hyppreferencelist)=cosh then hyp_preference:=cosh; let trig_normalize!*; p:=f; if keepalltriglist={} or directionlist={combine} or directionlist={compact} then <<let trig_standardize!*; p:=p+0;clearrules(trig_standardize!*);>>; if modelist neq {} then begin if first(modelist)=trig then <<let hyp2trig!*; p:=p+0;clearrules(hyp2trig!*);p:=behandle(p);>>; if first(modelist)=hyp then <<p:=behandle(p);let trig2hyp!*;p:=p+0;clearrules(trig2hyp!*)>>; if first(modelist)=expon then <<let trig2exp!*;p:=p+0;clearrules(trig2exp!*);>>; end; if first(directionlist)=expand then <<let trig_expand!*;p:=p+0;clearrules(trig_expand!*);>>; if first(directionlist)=combine then <<let trig_combine!*;p:=p+0;clearrules(trig_combine!*); if onlytan=1 and (keepalltriglist={keepalltrig}) then <<let subtan!*;p:=p+0;clearrules(subtan!*)>>;>>; clearrules(trig_normalize!*); if first(directionlist)=compact then begin % load compact; % Loaded at beginning. let trig_expand!*;p:=p+0;clearrules(trig_expand!*); p:=compact(f,{sin(x)**2+cos(x)**2=1}); end; end; return p; end; procedure degree(p,x); begin scalar h; if p=0 then h:=inf else h:=deg(num(p),x)-deg(den(p),x); return h; end; procedure balanced(p,x); if deg(num(p),x)=2*deg(den(p),x) then 1 else 0; procedure coordinated(p,x); begin scalar k,pneu,e,o,j; k:={}; e:=0; o:=0; pneu:=num(p); for j:=0:deg(pneu,x) do <<if coeffn(pneu,x,j) neq 0 then k:=j.k>>; for j:=1:length(k) do <<if fixp(part(k,j)/2) then e:=1 else o:=1>>; if o=e then return 0 else return 1; end; procedure trig2ord(p,x,y); begin if balanced(p,x) neq 1 or balanced(p,y) neq 1 then write "error using trig2ord: polynomial not balanced."; if coordinated(p,x) neq 1 or coordinated(p,y) neq 1 then write "error using trig2ord: polynomial not coordinated"; return sub(x=sqrt(x),y=sqrt(y),x**degree(p,x)*y**degree(p,y)*p); end; procedure ord2trig(p,x,y); x**(-degree(p,x))*y**(-degree(p,y))*sub(x=x**2,y=y**2,p); procedure factor_trig_poly(p,x,y); begin scalar j,p1,flist1,flist,d; p1:=trig2ord(p,x,y); d:=den(p1); flist1:=factorize(num(p1)); flist:={}; for j:=1:length(flist1) do flist:=ord2trig(part(flist1,j),x,y).flist; if d neq 1 then flist:=(1/d).flist; return flist; end; procedure subpoly2trig(p,x); begin scalar r,d; d:=degree(den(p),x); r:=p*x**d; r:=sub(x=cos(x)+i*sin(x),r); r:=r*(cos(x)-i*sin(x))**d; return r; end; procedure subpoly2hyp(p,x); begin scalar r,d; d:=degree(den(p),x); r:=p*x**d; r:=sub(x=cosh(x)+sinh(x),r); r:=r*(cosh(x)-sinh(x))**d; return r; end; procedure varget(p); begin scalar q,l,h; q:=factorize(p); h:=0; for each l in q do <<if not(numberp(l)) then begin if (h=0) and length(l)=1 then h:=l else h:=1; end; >>; if h=0 then h:=1; return h; end; procedure numberget(p); begin scalar q,d,l,h; q:=factorize(p); d:=1; h:=0; for each l in q do if numberp(l) then d:=d*l; return d; end; procedure argumenttriglist(term1, list1); begin scalar head1,j,l; l:=arglength(term1); if l<1 then return list1; head1:=part(term1,0); if head1 = sinh or head1 = cosh or head1 = sin or head1 = cos then list1:=append(list1,{part(term1,1)}) else for j:=1:l do list1:=argumenttriglist(part(term1,j),list1); return list1; end; procedure triggcd(p,q,x); begin scalar p1,q1,g1,g,u,d,nu,h,complex!*,err,l; on complex; nu:=numberget(x); err:=0; if varget(x)=1 then err:=1 else begin l:=argumenttriglist(p,{}); for d:=1:length(l) do if not(fixp(df(part(l,d),varget(x))/nu)) and not(freeof(part(l,d),varget(x))) then err:=1; l:=argumenttriglist(q,{}); for d:=1:length(l) do if not(fixp(df(part(l,d),varget(x))/nu)) and not(freeof(part(l,d),varget(x))) then err:=1; end; if err=0 then begin p1:=trigsimp1(p,{}); p1:=sub(sin(varget(x))=sin(varget(x)/nu), cos(varget(x))=cos(varget(x)/nu),sinh(varget(x)) =sinh(varget(x)/nu), cosh(varget(x))=cosh(varget(x)/nu),p1); p1:=trigsimp1(p1,{}); q1:=trigsimp1(q,{}); q1:=sub(sin(varget(x))=sin(varget(x)/nu), cos(varget(x))=cos(varget(x)/nu),sinh(varget(x)) =sinh(varget(x)/nu), cosh(varget(x))=cosh(varget(x)/nu),q1); q1:=trigsimp1(q1,{}); p1:=sub(sin(varget(x))=(xx_x-1/xx_x)/(2i), cos(varget(x)) =>xx_x/2+1/(2xx_x), sinh(varget(x))=(yy_y-1/yy_y)/2, cosh(varget(x))=yy_y/2+1/(2yy_y),p1); q1:=sub(sin(varget(x))=(xx_x-1/xx_x)/(2i), cos(varget(x)) =>xx_x/2+1/(2xx_x), sinh(varget(x))=(yy_y-1/yy_y)/2, cosh(varget(x))=yy_y/2+1/(2yy_y),q1); if balanced(p1,xx_x)+balanced(q1,xx_x)+coordinated(p1,xx_x)+ coordinated(q2,xx_x)+balanced(p1,yy_y)+balanced(q1,yy_y)+ coordinated(p1,yy_y)+coordinated(q2,yy_y) neq 8 then d:=1 else begin p1:=trig2ord(p1,xx_x,yy_y); q1:=trig2ord(q1,xx_x,yy_y); g1:=gcd(num(p1),num(q1)); g:=ord2trig(g1,xx_x,yy_y)/lcm(den(p1),den(p2)); h:=subpoly2trig(g,xx_x); h:=subpoly2hyp(h,yy_y); h:=sub(xx_x=varget(x)*nu,yy_y=varget(x)*nu,h); h:=trigsimp1(h,{}); h:=factorize(num(h)); d:=1; for each r in h do if not(numberp(r)) then d:=d*r; end; end else d:="error using triggcd, basis not possible."; return d; end; procedure trigfactorize(p,x); begin scalar l,q,f,r,d,h,s,u,err,complex!*; on complex; nu:=numberget(x); err:=0; if varget(x)=1 then err:=1 else begin l:=argumenttriglist(p,{}); for d:=1:length(l) do if not(fixp(df(part(l,d),varget(x))/nu)) and not(freeof(part(l,d),varget(x))) then err:=1; end; if err=1 then rederr("error using trigfactorize, basis not possible") else begin q:=trigsimp1(p,{}); q:=sub(sin(varget(x))=sin(varget(x)/nu), cos(varget(x))=cos(varget(x)/nu),sinh(varget(x)) =sinh(varget(x)/nu), cosh(varget(x))=cosh(varget(x)/nu),q); q:=trigsimp1(q,{}); q:=sub(sin(varget(x))=(xx_x-1/xx_x)/(2i), cos(varget(x)) =xx_x/2+1/(2xx_x), sinh(varget(x))=(yy_y-1/yy_y)/2,cosh(varget(x)) =yy_y/2+1/(2yy_y),q); if balanced(q,xx_x)+coordinated(q,xx_x)+ balanced(q,yy_y)+coordinated(q,yy_y)<4 then f:={p} else begin q:=factor_trig_poly(q,xx_x,yy_y); f:={}; d:=1; for each r in q do <<h:=subpoly2trig(r,xx_x); h:=subpoly2hyp(h,yy_y); h:=sub(xx_x=varget(x)*nu,yy_y=varget(x)*nu,h); h:=trigsimp1(h,{}); if freeof(h,varget(x)) then d:=d*h else begin for each u in factorize(h) do <<if freeof(u,varget(x)) then begin d:=d*u; h:=h/u; end; >>; f:=reverse(h.reverse(f)); end; >>; if d neq 1 then f:=d.f; end; end; return f; end; symbolic procedure trigsimp!*(f); trigsimp1(reval car f,'list.for each w in cdr f collect reval w); symbolic put('trigsimp,'psopfn, 'trigsimp!*); endmodule; end;