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