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; 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