%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Neil Langmead November 1996 ZIB Berlin
% routines to evaluate trigonometric integrals
%
% main routine is trigint
% substitution variable is always u
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
module trigint;
create!-package ('(trigint),nil);
global '(!*tracetrig);
switch tracetrig;
off tracetrig; % off by default;
algebraic;
load_package limits;
on usetaylor;
load_package misc;
% need a knowledge base of the substitutions
expr procedure sub_a(exp,var);
sub({sin(var)=2*u/(1+u^2), cos(var)=(1-u^2)/(1+u^2)},exp);
expr procedure sub_b(exp,var);
sub({sin(var)=(u^2-1)/(u^2+1),
cos(var)=2*u/(u^2+1)},exp);
expr procedure sub_c(exp,var);
sub({sin(var)=2*u/(1+u^2), cos(var)=(u^2-1)/(1+u^2)},exp);
expr procedure sub_d(exp,var);
sub({sin(var)=u/(sqrt(1+u^2)), cos(var)=1/(sqrt(1+u^2))},exp);
% applying the substitutions to their integrals
expr procedure apply_a(exp,var);
begin scalar answer, result;
answer:=sub_a(exp,var);
answer:=answer* 2/(1+u^2);
result:=int(answer,u);
result:=sub({u=tan(var/2)},result);
result:=result+k(result,var,pi)*floor((var-pi)/(2*pi));
return result;
end;
expr procedure apply_b(exp,var);
begin scalar answer, result;
answer:=sub_b(exp,var); answer:=answer*2/(1+u^2);
result:=int(answer,u); result:= sub({u=tan(var/2+pi/4)}, result);
result:=result+k(result,var,pi/2)*floor((var-pi/2)/(2*pi));
return result;
end;
expr procedure apply_c(exp, var);
begin scalar answer, result;
answer:=sub_c(exp,var); answer:= answer*(-2/(1+u^2));
result:=int(answer,u); result:=sub({u=1/(tan(var/2))},result);
result:=result+k(result,var,0)*floor(var/(2*pi));
return result;
end;
expr procedure apply_d(exp, var);
begin scalar answer, result;
answer:=sub_d(exp,var);
answer:=answer*(1/(1+u^2));
result:=int(answer,u); result:= sub({u=tan(var)},result);
result:=result+k(result,var,pi/2)*floor((var-pi/2)/pi);
return result;
end;
% some trigonometric substitutions which occur very frequently
trig_rules:= {
(sec(~x))^2 => 1/(cos(x))^2,
%tan(~x) => sin(x)/cos(x),
(tan(~x))^2 => (sec(x))^2-1
%2*sin(~x/2)*cos(~x/2) => sin(x)
%(cos(~x))^2 => 1-(sin(x))^2 % (sin(~x))^2 => 1-(cos(x))^2
};
%let trig_rules;
%for all x let sin(x)^2+cos(x)^2=1,
% 2*sin(x/2)*cos(x/2)=sin(x);
% procedure to see if integral is returned unevaluated
expr procedure unevalp(exp);
begin scalar finished, k; k:=0; finished:=0;
while ( finished=0 and k<=arglength(exp)) do
<<
if(part(exp,k)=int) then finished:=1 else k:=k+1;
>>;
if (finished=1) then return t else nil;
end;
% procedure to evaluate K
expr procedure k(exp,var,val);
limit!-(exp,var,val)-limit!+(exp,var,val);
% two routines to see if we have either unevaluated limits or ints in our
% result. If we have, then try a different substitution
expr procedure uneval_int(exp);
begin scalar temp;
if( not freeof(exp,int)) then return temp:=t else temp:=nil;
return temp;
end;
expr procedure uneval_lim(exp);
begin scalar temp;
if(not freeof(exp,limit!-)) then return t else
<<
if(not freeof(exp,limit!+)) then return t else nil;
>>;
end;
expr procedure fail_a(exp,var);
begin scalar temp;
temp:=apply_a(exp,var);
if(uneval_lim(temp)) then return t else
<<
if(uneval_int(temp)) then return t
else return nil;
>>;
end;
expr procedure fail_b(exp,var);
begin scalar temp;
temp:=apply_b(exp,var);
%temp:=temp+k(temp,var,pi/2)*floor((var-pi/2)/(2*pi));
if(uneval_lim(temp)) then return t else
<<
if(uneval_int(temp)) then return t else return nil;
>>;
end;
expr procedure fail_c(exp,var);
begin scalar temp;
temp:=apply_c(exp,var);
%temp:=temp+k(temp,var,0)*floor(var/(pi));
if(uneval_lim(temp)) then return t else
<<
if(uneval_int(temp)) then return t else return nil;
>>;
end;
expr procedure fail_d(exp,var);
begin scalar temp;
temp:=apply_d(exp,var);
%temp:=temp+k(temp,var,pi/2)*floor((var-pi/2)/pi);
if(uneval_lim(temp)) then return t else
<<
if(uneval_int(temp)) then return t else return nil;
>>;
end;
expr procedure fail(exp);
if(uneval_lim(exp)) then t else
<< if(uneval_int(exp)) then t else nil; >>;
let log(-1) => i*pi; % really important. If further log rules are needed,
% take a look at the ratint package, and the module
% convert, which contains an extensive list of such
% rules
expr procedure trigint(exp,var);
begin scalar answer, answer_1, answer_2, answer_3, answer_4, result;
%off mesgs; % off by default
% check for correct input
if(freeof(exp,sin) and freeof(exp,cos) and freeof(exp,tan)) then <<
if(lisp !*tracetrig) then
write "expression free of sin, cos tan, proceeding with standard integration";
return int(exp,x); >>;
on usetaylor;
if freeof(exp,sin(var)) then % we use substitution (a)
<<
answer:=apply_a(exp,var);
%answer:=answer+k(answer,var,pi)*floor((var-pi)/(2*pi));
if(fail(answer)) then % system can't evaluate after subs
<<
if(lisp !*tracetrig) then
write "system can't integrate after substitution A,
trying again";
answer_2:=apply_b(exp,var);
if(fail(answer_2)) then
<<
if(lisp !*tracetrig) then write "trying again with substitution B";
answer_3:=apply_c(exp,var);
if(fail(answer_3)) then
<<
if(lisp !*tracetrig) then
write "and again with substitution C";
answer_4:=apply_d(exp,var);
if(fail(answer_4)) then
<<
if(lisp !*tracetrig) then
write "failed in all attempts, system cannot integrate";
return answer;
>> else return answer_4;
>> else return answer_3;
>> else return answer_2;
>> else return answer;
%let trig_rules;
%if(unevalp(answer)) then rederr "system cannot integrate after subs"
%else nil;
>>
else
% we use substitution b,c or d
<<
if(freeof(exp,cos(var))) then % use substitution b
<<
answer:=apply_b(exp,var);
if(fail(answer)) then
<<
if(lisp !*tracetrig) then write "failed with substitution B: system could not
integrate after subs, trying A";
answer_2:=apply_a(exp,var);
if(fail(answer_2)) then
<<
if(lisp !*tracetrig) then write "failed with A: trying C now";
answer_3:=apply_c(exp,var);
if(fail(answer_3)) then
<<
if(lisp !*tracetrig) then write "failed with C: trying D now";
answer_4:=apply_d(exp,var);
if(lisp !*tracetrig) then
write "trying all possible substitutions";
if(fail(answer_4)) then rederr "system can't integrate after
subs"
else return answer_4;
>> else return answer_3;
>> else return answer_2;
>> else return answer;
>>
else <<
% now describe situations best for (c) and (d) G and R sect 2.504
if(sub({sin(var)=-sin(var),cos(var)=-cos(var)},exp)=exp) then
% d is the best sub in this case
<<
if(lisp !*tracetrig) then write
"using heuristics: G & R section 2.504 to integrate ";
answer:=apply_d(exp,var);
if(fail(answer)) then
<<
if(lisp !*tracetrig) then write "subs D failed, trying now with A";
answer_2:=apply_a(exp,var);
if(fail(answer_2)) then
<<
if(lisp !*tracetrig) then write "subs B falied, trying with sub C";
answer_3:=apply_b(exp,var);
if(fail(answer_3)) then
<<
if(lisp !*tracetrig) then write "sub C falied, trying sub D";
answer_4:=apply_c(exp,var);
if(fail(answer_4)) then rederr "can't integrate after subs"
else return answer_4;
>> else return answer_3;
>> else return answer_2;
>> else return answer;
>>
else << % no guidelines, try each substitution in turn, and return the
% best possible answer, if there is one
answer:=apply_a(exp,var);
if(fail(answer)) then
<<
if(lisp !*tracetrig) then write "not using heuristics,
attempting subs in order: trying A";
answer_2:=apply_b(exp,var);
if(fail(answer_2)) then
<<
if(lisp !*tracetrig) then write "A failed, trying B";
answer_3:=apply_c(exp,var);
if(fail(answer_3)) then
<< answer_4:=apply_d(exp,var);
if(fail(answer_4)) then rederr "can't do it"
else return answer_4;
>> else return answer_3;
>> else return answer_2;
>> else return answer;
>>;
>>; % for the else just before the comment on G and R
>>;
return answer;
end;
endmodule;
end;