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(a)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))) /
> 2)
> + 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