To be honest, the revised Complex_SqRt function you present is a bit too complicated for me. <G>
So, my thinking is that, instead of re-writing the algorithm of Complex_SqRt and similar routines, it would be simpler and less error-prone in these situations, to utilize a new, alternative SqRt function, eg. :
Function Clipped_SqRt (x : double) : double; {Alternative SqRt : Clips negative arguments, returning zero} begin {Clipped_SqRt} if x > 0 then Clipped_SqRt := SqRt (x) else Clipped_SqRt := 0 end; {Clipped_SqRt}
Of course, you would only apply this alternative to SqRt when the argument cannot theoretically be negative ...
Joe.
-----Original Message----- From: Emil Jerabek [SMTP:ejer5183@artax.karlin.mff.cuni.cz] Sent: Wednesday, January 09, 2002 3:54 AM To: gpc@gnu.de Subject: Complex_SqRt bug
The complex square root function crashes with a funny error message on some arguments close to the real line:
[pas]% cat sqrt.pas program SqrtTest (Output);
var C: Complex;
begin C := Sqrt (Cmplx (4.2, 1e-8)); WriteLn ('OK') end. [pas]% uname -a ; gpc -v Linux localhost.localdomain 2.4.2-2 #1 Sun Apr 8 19:37:14 EDT 2001 i586 unknown Reading specs from /usr/local/lib/gcc-lib/i586-pc-linux-gnu/2.95.2/specs gpc version 20011222, based on 2.95.2 19991024 (release) [pas]% gpc sqrt.pas [pas]% ./a.out ./a.out: argument to `SqRt' is < 0 (error #708 at 805d837) [pas]%
The square root is computed as
\sqrt{X+iY} = \sqrt{(A+X)/2} + i \sqrt{(A-X)/2} sgn(Y),
where A=|X+iY|. It seems that due to inaccuracies etc., A may turn out to be _less_ than |X|, if |Y| << |X|, which means that SqRt((A+-X)/2) gets a negative argument. I suggest the following patch.
diff -U3 orig/p/rts/math.pas p/rts/math.pas --- orig/p/rts/math.pas Mon Sep 3 12:48:15 2001 +++ p/rts/math.pas Mon Jan 7 21:34:43 2002 @@ -192,8 +192,7 @@
function Complex_SqRt (z: Complex): Complex; var
- a: Double;
- Sign: Integer;
- a, b: Double;
begin { SqRt (x + yi) = +/- (SqRt ((a + x) / 2) +/- i * SqRt ((a - x) / 2))
@@ -203,10 +202,15 @@ Principal value defined in the Extended Pascal standard: Exp (0.5 * Ln (z)), i.e. Sign (Re (SqRt (z)) := 1, Sign (Im (SqRt (z))) := Sign (Im (z)) } a := Abs (z);
- if Im (z) < 0 then Sign := -1 else Sign := 1;
- if a <> 0
- then Complex_SqRt := Cmplx (SqRt ((a + Re (z)) / 2), Sign * SqRt ((a
- Re (z)) / 2))
- else Complex_SqRt := 0
- if a = 0 then Complex_SqRt := 0
- else begin
- b := SqRt ((a + Abs (Re (z))) / 2);
- if Abs (Im (z)) > Abs (Re (z)) then a := SqRt ((a - Abs (Re (z))) /
else a := Abs (Im (z)) / (2 * b);- if Re (z) < 0 then begin var c: Double = a; a := b; b := c end;
- if Im (z) < 0 then Complex_SqRt := Cmplx (b, -a)
else Complex_SqRt := Cmplx (b, a)- end
end;
function Complex_Ln (z: Complex): Complex;
It would also work to simply insert
if Abs (Re (z)) > a then a := Abs (Re (z));
into the original version, but the patch above gives more accurate results (even for arguments which do not trigger the bug).
When talking about math.pas, I also suggest a little optimization to the complex logarithm function; \log \sqrt x = (\log x)/2, but the RHS is (I guess) computed faster and more precisely:
diff -U3 orig/p/rts/math.pas p/rts/math.pas --- orig/p/rts/math.pas Mon Sep 3 12:48:15 2001 +++ p/rts/math.pas Mon Jan 7 22:10:08 2002 @@ -211,7 +211,8 @@
function Complex_Ln (z: Complex): Complex; begin
- Complex_Ln := Cmplx (Ln (Abs (z)), ArcTan2 (Im (z), Re (z)))
- Complex_Ln := Cmplx (Ln (Sqr (Re (z)) + Sqr (Im (z))) /2,
ArcTan2 (Im (z), Re (z)))end;
function Complex_Exp (z: Complex): Complex;
Emil Jerabek