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
Emil Jerabek wrote:
The following patch to p/rts/math.pas solves these problems (I hope it doesn't introduce others).
Before I look closer into the routine (I'm quite short of time at the moment), one remark:
function Log1P (x: Real): Real; asmname 'log1p';
According to the manpage:
: CONFORMING TO : BSD
So this function might not be available everywhere, and we need a replacement for other systems. Can you please provide one (written in Pascal)? I can deal with autoconf then (i.e., checking if the function is available and using the replacement if not).
Frank
Frank wrote:
function Log1P (x: Real): Real; asmname 'log1p';
According to the manpage:
: CONFORMING TO : BSD
Oops. My manpage says the same. OTOH, info libc doesn't say anything (it usually does for non-ISO C89 functions), the header says
#if defined __USE_MISC || defined __USE_XOPEN_EXTENDED || defined __USE_ISOC9X
("MISC" means the common intersection of BSD and System V), and I've succesfully compiled a simple test program even with -D__STRICT_ANSI__. This seems like quite a mess to me.
So this function might not be available everywhere, and we need a replacement for other systems. Can you please provide one (written in Pascal)? I can deal with autoconf then (i.e., checking if the function is available and using the replacement if not).
Here it is. The C library function is certainly preferable when present -- it is likely to have a hardware support, so it is faster (although I tried to optimize the Pascal implementation for speed, not code clarity).
--------------------------------------- var ExpEpsReal: Integer; Log1PCoefs: array [0 .. BitSizeOf (Real) div 3] of Real; DummyMantissa: LongReal;
function Log1P (X: Real) = Y: Real; var I, N: Integer; Z, U: Real; begin if (X <= -0.5) or (X >= 1) then Y := Ln (1 + X) else begin U := X / (2 + X); Z := Sqr (U); if Z <= MinReal then Y := X else begin SplitReal (Z, N, DummyMantissa); N := ExpEpsReal div N; Y := 0; for I := N downto 0 do Y := Y * Z + Log1PCoefs[I]; Y := Y * U end end end;
var I: Integer;
to begin do begin SplitReal (EpsReal, ExpEpsReal, DummyMantissa); Dec (ExpEpsReal); for I := 0 to BitSizeOf (Real) div 3 do Log1PCoefs[I] := 2 / (I * 2 + 1) end; -----------------------------------------------
Before I look closer into the routine (I'm quite short of time at the moment), one remark:
To save your time, here is some mathematics behind my Complex_Arctan reimplementation:
Let z = x+iy. Arctan (z) is defined as -i/2 Ln((1+iz)/(1-iz)). Thus
Re(arctan z) = 1/2 arg ((1+iz)/(1-iz)) = = 1/2 arg ( (1+iz)(1+iz') / |1-iz|^2 ) = = 1/2 arg ((1+iz)(1+iz')) = = 1/2 arg (1+i(z+z')-|z|^2) = = 1/2 arg (1-x^2-y^2+2ix) = = 1/2 Arctan2 (2x, (1+y)(1-y)-x^2)
where z' denotes the conjugate of z, and
Im(arctan z) = -1/2 ln |(1+iz)/(1-iz)| = = -1/2 ln |1+iz|/|1-iz| = = 1/2 ln |1-iz|/|1+iz| = = 1/2 ln \sqrt{((1+y)^2 + x^2) / ((1-y)^2 + x^2)} = = 1/4 ln (((1+y)^2 + x^2) / ((1-y)^2 + x^2)) = = 1/4 ln (1 + 4y / ((1-y)^2 + x^2))
The function deals with |y| instead of y to save some "if"s; it's clear from the equations above that arctan(z') = (arctan z)'.
The "if c <= MinReal" branch (needed to avoid some overflows) is executed when |y|=1 and x^2 <= MinReal. In this case,
1/4 ln (((1+|y|)^2 + x^2) / ((1-|y|)^2 + x^2) = 1/4 ln ((4+x^2) / x^2) = [not exactly, but recall that x^2 << 1] 1/4 ln (4 / x^2) = 1/2 ln (2/x)
The Pascal version of log1p() uses the power series
ln ((1+u) / (1-u)) = 2 \sum_{n=0}^{\infty} u^{2n+1} / (2n+1) (|u|<1)
together with ln (1+x) = ln ((1+u) / (1-u)), where u = x/(2+x).
Emil Jerabek
Emil Jerabek wrote:
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))
With the improved Abs (and therefore Ln), this doesn't seem to happen anymore. But I wrote a test prorgam (emil9.pas) which shows some cases of problems even with the improved Abs. They all work ok with your code. I'm attaching the test program -- you might want to take a look at it and possible improve or extend it.
The following patch to p/rts/math.pas solves these problems (I hope it doesn't introduce others).
Thanks for the code. I'm putting it in the RTS.
To save your time, here is some mathematics behind my Complex_Arctan reimplementation:
[...]
Mathematically it all seems alright. I didn't try to prove the numeric stability, but my tests show that it's better than that of the previous version. Just one thing:
The "if c <= MinReal" branch (needed to avoid some overflows) is executed when |y|=1 and x^2 <= MinReal. In this case,
From the code, the condition for this case is:
| 1 - |y| | <= |x| and x^2 <= MinReal
To conclude |y| = 1, you need that EpsReal^2 >= MinReal, don't you? While this might be true on every FP system, would it cause problems in your code if not?
Frank
Frank wrote:
Mathematically it all seems alright. I didn't try to prove the numeric stability, but my tests show that it's better than that of the previous version. Just one thing:
The "if c <= MinReal" branch (needed to avoid some overflows) is executed when |y|=1 and x^2 <= MinReal. In this case,
From the code, the condition for this case is:
| 1 - |y| | <= |x| and x^2 <= MinReal
To conclude |y| = 1, you need that EpsReal^2 >= MinReal, don't you? While this might be true on every FP system, would it cause problems in your code if not?
Frank
The "c <= MinReal" test is completely wrong :-( It should be "c <= 5 / MaxReal".
But back to your question. The code really assumes EpsReal^2 > 5/MaxReal. If you feel that this may be a problem, something more complicated is needed:
if c + Sqr (b) <= 5 / MaxReal then b := -Ln (Abs (Cmplx (Re (z), b)) / (1 + Abs (Im (z)))) / 2 {*} else if Abs (b) > Abs (Re (z)) then b := Log1P (4 * (Abs (Im (z)) / b) / (b + Re (z) * (Re (z) / b))) / 4 else b := Log1P (4 * (Abs (Im (z)) / Re (z)) / (b * (b / Re (z)) + Re (z))) / 4;
This still assumes 1/EpsReal < 1.6*MaxReal. I can't imagine a floating-point implementation violating this; anyway, in such a case, the line with {*} above should be replaced by
b := -Ln (Abs (Cmplx (Re (z), b)) / Abs (Cmplx (Re (z), 1 + Abs (Im (z))))) / 2
This code is essentially untested (i.e. it works on i586/Linux, but that's not a target with EpsReal^2 <= 5/MaxReal).
Emil Jerabek