module result;
% Author: James H. Davenport.
fluid '( !*rationalize !*tra gaussiani !*trmin intvar );
exports combine!-logs;
symbolic procedure combine!-logs(coef,logarg);
% Attempts to produce a "simple" form for COEF*(LOGARG). COEF is
% prefix, LOGARG an SQ (and already log'ged for historical reasons).
begin
scalar ans,dencoef,parts,logs,lparts,!*rationalize,trueimag;
!*rationalize:=t; % A first attempt to use this technology.
coef:=simp!* coef;
parts:=split!-real!-imag numr coef;
if null numr cdr parts then return multsq(coef,logarg);
% Integrand was, or seemed to be, purely real.
dencoef:=multf(denr coef,denr logarg);
if !*tra then <<
printc "attempting to find 'real' form for";
mathprint list('times,list('plus,prepsq car parts,
list('times,prepsq cdr parts,'i)),
prepsq logarg) >>;
logarg:=numr logarg;
logs:= 1 ./ 1;
while pairp logarg do <<
if ldeg logarg neq 1 then interr "what a log";
if atom mvar logarg then interr "what a log";
if car mvar logarg neq 'log then interr "what a log";
logs:=!*multsq(logs,
!*exptsq(simp!* argof mvar logarg,lc logarg));
logarg:=red logarg >>;
logs:=rationalizesq logs;
ans:=multsq(!*multsq(car parts,logs),1 ./ dencoef); % real part
% Now to apply theory i*log(a+i*b) = atan(a/b) + (i/2 log (a^2+b^2))
lparts:=split!-real!-imag numr logs;
if numr difff(denr cdr lparts,intvar)
then interr "unexpected denominator";
lparts:=!*multsq(denr cdr lparts ./ 1,car lparts) . cdr lparts;
if not onep denr car lparts then interr "unexpected denominator";
% We have discarded the logarithm of a constant: good riddance
trueimag:=quotsq(addf(!*exptf(numr car lparts,2),
!*exptf(numr cdr lparts,2)) ./ 1,
!*exptf(denr logs,2) ./ 1);
if numr diffsq(trueimag,intvar)
then ans:=!*addsq(ans,
!*multsq(gaussiani ./ multf(2,dencoef),
!*multsq(simplogsq trueimag,cdr parts)));
trueimag:=!*multsq(car lparts,!*invsq(numr cdr lparts ./ 1));
if numr diffsq(trueimag,intvar)
then ans:=!*addsq(ans,!*multsq(!*multsq(cdr parts,1 ./ dencoef),
!*k2q list('atan,prepsq!* trueimag)));
return ans;
end;
symbolic procedure split!-real!-imag sf;
% Returns coef real.imag as SQs.
if null sf then (lambda z; z . z) (nil ./ 1)
else if numberp sf then (sf ./ 1) . (nil ./ 1)
else if domainp sf then interr "can't handle arbitrary domains"
else begin
scalar cparts,rparts,mv,tmp;
cparts:=split!-real!-imag lc sf;
rparts:=split!-real!-imag red sf;
mv:=split!-real!-imagvar mvar sf;
if null numr cdr mv % main variable totally real
then <<
tmp:= lpow sf .* 1 .+ nil ./ 1;
return !*addsq(!*multsq(car cparts,tmp),car rparts) .
!*addsq(!*multsq(cdr cparts,tmp),cdr rparts) >>;
if null numr car mv then <<
mv:=!*exptsq(cdr mv,ldeg sf);
% deal with powers of i
if not evenp(ldeg sf / 2) then mv:=negsq mv;
if not evenp ldeg sf
then return !*addsq(!*multsq(negsq cdr cparts,mv),car rparts) .
!*addsq(!*multsq(car cparts,mv),cdr rparts)
else return !*addsq(!*multsq(car cparts,mv),car rparts) .
!*addsq(!*multsq(cdr cparts,mv),cdr rparts) >>;
% Now we have to handle the general case.
cparts:=mul!-complex(cparts,exp!-complex(mv,ldeg sf));
return !*addsq(car cparts,car rparts) .
!*addsq(cdr cparts, cdr rparts)
end;
symbolic procedure mul!-complex(a,b);
!*addsq(!*multsq(negsq cdr a,cdr b),!*multsq(car a,car b)) .
!*addsq(!*multsq(car a,cdr b),!*multsq(cdr a,car b));
symbolic procedure exp!-complex(a,n);
if n=1 then a
else if evenp n then exp!-complex(mul!-complex(a,a),n/2)
else mul!-complex(a,exp!-complex(mul!-complex(a,a),n/2));
symbolic procedure split!-real!-imagvar mv;
% Returns a pair of sf.
if mv eq 'i then (nil ./ 1) . (1 ./ 1)
else if atom mv then (mv .** 1 .* 1 .+ nil ./ 1) . (nil ./ 1)
else if car mv eq 'sqrt then begin
scalar n,rp,innersqrt,c;
n:=simp!* argof mv;
if denr n neq 1 then interr "unexpected denominator";
rp:=split!-real!-imag numr n;
if null numr cdr rp and minusf numr car rp and
null involvesf(numr car rp,intvar) then
return (nil ./ 1) . simpsqrtsq negsq car rp;
if null numr cdr rp
then return (mv .** 1 .* 1 .+ nil ./ 1) . (nil ./ 1);
% totally real.
% OK - it's a general (ish) complex number A+iB
% its square root is going to be C+iD where
% C^2 = (A+sqrt(A^2+B^2))/2 (+ve sign of sqrt to make C +ve)
% C is square root of this
% D is C * (sqrt(A(2+B^2) -A)/B
% Note that D has a non-trivial denominator. We could avoid this at
% the cost of generating non-independent square roots (yuck).
% Note that the above checks have ensured this den. is non-zero.
if numr car rp then
innersqrt:=simpsqrtsq !*addsq(!*exptsq(car rp,2),
!*exptsq(cdr rp,2))
else innersqrt:=cdr rp; % pure imaginary case
c:=simpsqrtsq multsq(!*addsq(car rp, innersqrt), 1 ./ 2);
return c . !*multsq(!*multsq(c,!*invsq cdr rp),
!*addsq(innersqrt,negsq car rp));
end
else (mv .** 1 .* 1 .+ nil ./ 1) . (nil ./ 1);
% What the heck: pretend it's real.
endmodule;
end;