Artifact 7f11a614fc40d13e85fded9d1abf981373eee19da6a69f84735af817f016105c:
- Executable file
r37/packages/algint/log2atan.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: 5741) [annotate] [blame] [check-ins using] [more...]
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;