Please see below ...
Joe.
-----Original Message----- From: Emil Jerabek [SMTP:ejer5183@artax.karlin.mff.cuni.cz] Sent: Friday, January 25, 2002 5:23 AM To: gpc@gnu.de Subject: Complex_Arctan
Hi all,
I've had a look at the complex arctangent function in the RTS. It turned out that
- it raises a runtime error or returns an infinite imaginary part on
arguments too close to the singularities +/- i, for example ArcTan (Cmplx (1e6 * MinReal, +/- 1))
- it often returns a real number when it should return something with a
small imaginary part instead (example: ArcTan (Cmplx (1e17, 1e17)) on IA32, YMMV), and in general, the imaginary part of the result is very unreliable when |Im (z)| << |Re (z)|^2 (because it is essentially computed as `Ln (something close to 1)'; example: Im (Arctan (Cmplx (1000, 0.0001))) gives 9.999989726e-11 instead of 9.999990000e-11, again this is machine-dependent)
The following patch to p/rts/math.pas solves these problems (I hope it doesn't introduce others).
function Log1P (x: Real): Real; asmname 'log1p';
function Complex_Arctan (z: Complex): Complex; const Ln2 = 0.6931471805599453094172321214581765680755; var b, c, d: Real; begin c := Sqr (Re (z)); b := 1 - Abs (Im (z)); d := ((1 + Abs (Im (z))) * b - c) / 2;
[Joe da Silva]
I wonder, would it be better to use the following instead, to reduce problems with overflow for large values (magnitude) of Re(z) and/or Im(z)? :
c := Re(z) * 0.5 * Re(z); b := 1 - Abs(Im(z)); d := (1 + Abs(Im(z))) * 0.5 * b - c;
Also, as a general question, is division more expensive than multiplication, or does the complier's optimization take care of such matters anyway?
if Abs (b) > Abs (Re (z)) then b := Log1P (4 * (Abs (Im (z)) / b) / (b + Re (z) * (Re (z) / b))) / 4 else if c <= MinReal then b := (Ln2 - Ln (Abs (Re (z)))) / 2 else b := Log1P (4 * (Abs (Im (z)) / Re (z)) / (b * (b / Re (z)) + Re (z))) / 4; if Im (z) < 0 then b := -b; Complex_Arctan := Cmplx (ArcTan2 (Re (z), d) / 2, b) end;
Emil Jerabek