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(a)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
>