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