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
Emil Jerabek wrote:
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.
Thanks (emil3.pas). Though I can't reproduce the problem, even though I'm also on IA32, your patch looks good. (I vaguely recall hearing about numerically stable algorithms once ...)
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:
Also applied, thanks. If you find other things to improve in the maths algorithms, just go ahead (numerics has never been my favourite area ;-).
Frank
Frank Heckenbach wrote:
Emil Jerabek wrote: If you find other things to improve in the maths algorithms, just go ahead (numerics has never been my favourite area ;-).
While reading this thread, I had a look to math.pas and found an other well known issue. The function
Complex_Abs := SqRt (Sqr (Re (z)) + Sqr (Im (z)))
overflows if either Re or Im is in between sqrt(MaxReal) (~10^154) and MaxReal (~10^308), while the result is (usually) below MaxReal. One Quick and Dirty trick is adapted from the Numerical Recipes library (version 2), where the function is called pythag:
function Complex_Abs (z: Complex): Double; var absx, absy: double; begin absx:=abs(Re(z)); absy:=abs(Im(z)); if absx>absy then Complex_Abs:=absx*sqrt(1.0+sqr(absy/absx)) else if absy=0 then Complex_Abs:=0.0 else Complex_Abs:=absy*sqrt(1.0+sqr(absx/absy)); end;
Now I had a look for this problem, both in the djgpp libc and libm libraries, where there is a function hyphot which does the job, and is said to be ieee754 standard, and in the CEPHES long double math library, (which defines cabs), which I use (with an import unit) when I need accuracy, since djgpp has no long double library (only elementary + - / * operations I think).
libc overcomes the problem by using the extended accuracy and exponent range of the x87 floating point unit. gpc cannot use directly the same trick, since sqrt is computed in double accuracy only because djgpp does not define sqrtl (this is tested by configure). One way out would be to use an other configure test for hypot, to use directly this libc function (if this is possible).
libm and CEPHES use basically the same trick to avoid the need for extended exponent range: they extract the exponents, rescale the two arguments (by shifting their exponents) to be in acceptable range, compute and then rescale the result the opposite way. There are also auxiliary issues to ensure lowest bit accuracy described in %DJGPP%\src\libm\math\e_hypot.c (extracted from v2/djlsr203.zip). The libm hypot function could be also retrieved through configure test I think.
Basically what to do depends on what kind of complexity you are willing to incorporate into the compiler itself, as opposed to specialized libraries.
Maurice
Maurice Lombardi wrote:
Frank Heckenbach wrote:
Emil Jerabek wrote: If you find other things to improve in the maths algorithms, just go ahead (numerics has never been my favourite area ;-).
While reading this thread, I had a look to math.pas and found an other well known issue. The function
Complex_Abs := SqRt (Sqr (Re (z)) + Sqr (Im (z)))
overflows if either Re or Im is in between sqrt(MaxReal) (~10^154) and MaxReal (~10^308), while the result is (usually) below MaxReal.
Why not apply a little trigonometry.
Theta = arctan(Re/Im); C_Abs = Re * cos(theta);
with a few conditionals to avoid zeroes. Actually, the presence of zeroes greatly simplifies the computation.
You might reserve it for when either Re or Im is greater than sqrt(maxreal)/2 if you have qualms about the trigonometric function accuracy. Note the /2, so the sum doesn't overflow either.
(* only overflows if the result is an overflow value *) IF re = 0.0 THEN cabs := abs(im) ELSE IF im = 0.0 THEN cabs := abs(re) ELSE IF (abs(re) < crit) AND (abs(im < crit) THEN BEGIN (* sqrt of sum of squares, can't overflow *) END ELSE BEGIN (* trig version, can overflow, probably less accurate *) END;
CBFalconer wrote:
Maurice Lombardi wrote:
Frank Heckenbach wrote:
Emil Jerabek wrote: If you find other things to improve in the maths algorithms, just go ahead (numerics has never been my favourite area ;-).
While reading this thread, I had a look to math.pas and found an other well known issue. The function
Complex_Abs := SqRt (Sqr (Re (z)) + Sqr (Im (z)))
overflows if either Re or Im is in between sqrt(MaxReal) (~10^154) and MaxReal (~10^308), while the result is (usually) below MaxReal.
Why not apply a little trigonometry.
Theta = arctan(Re/Im); C_Abs = Re * cos(theta);
You mean /, not *?
with a few conditionals to avoid zeroes.
Not only zeros. Re/Im might overflow if Im is very small compared to Re.
You might reserve it for when either Re or Im is greater than sqrt(maxreal)/2 if you have qualms about the trigonometric function accuracy.
... and efficiency. I think one SqRt is quite a bit cheaper than two trigonomic functions.
Frank
Frank Heckenbach wrote:
CBFalconer wrote:
Maurice Lombardi wrote:
... snip ...
While reading this thread, I had a look to math.pas and found an other well known issue. The function
Complex_Abs := SqRt (Sqr (Re (z)) + Sqr (Im (z)))
overflows if either Re or Im is in between sqrt(MaxReal) (~10^154) and MaxReal (~10^308), while the result is (usually) below MaxReal.
Why not apply a little trigonometry.
Theta = arctan(Re/Im); C_Abs = Re * cos(theta);You mean /, not *?
with a few conditionals to avoid zeroes.
Not only zeros. Re/Im might overflow if Im is very small compared to Re.
You might reserve it for when either Re or Im is greater than sqrt(maxreal)/2 if you have qualms about the trigonometric function accuracy.
... and efficiency. I think one SqRt is quite a bit cheaper than two trigonomic functions.
You didn't read it all. It only comes into play when there is potential overflow.
CBFalconer wrote:
Frank Heckenbach wrote:
CBFalconer wrote:
Maurice Lombardi wrote:
... snip ...
While reading this thread, I had a look to math.pas and found an other well known issue. The function
Complex_Abs := SqRt (Sqr (Re (z)) + Sqr (Im (z)))
overflows if either Re or Im is in between sqrt(MaxReal) (~10^154) and MaxReal (~10^308), while the result is (usually) below MaxReal.
Why not apply a little trigonometry.
Theta = arctan(Re/Im); C_Abs = Re * cos(theta);You mean /, not *?
with a few conditionals to avoid zeroes.
Not only zeros. Re/Im might overflow if Im is very small compared to Re.
You might reserve it for when either Re or Im is greater than sqrt(maxreal)/2 if you have qualms about the trigonometric function accuracy.
... and efficiency. I think one SqRt is quite a bit cheaper than two trigonomic functions.
You didn't read it all.
Whoa, you know better than me what I read?
It only comes into play when there is potential overflow.
But Maurice's (or Numerical Recipes's) solution avoids the overflow as well with only one SqRt and some primitive operations (+, *, /, >). That's what I was comparing with.
Frank
Maurice Lombardi wrote:
Frank Heckenbach wrote:
Emil Jerabek wrote: If you find other things to improve in the maths algorithms, just go ahead (numerics has never been my favourite area ;-).
While reading this thread, I had a look to math.pas and found an other well known issue. The function
Complex_Abs := SqRt (Sqr (Re (z)) + Sqr (Im (z)))
overflows if either Re or Im is in between sqrt(MaxReal) (~10^154) and MaxReal (~10^308), while the result is (usually) below MaxReal. One Quick and Dirty trick is adapted from the Numerical Recipes library (version 2), where the function is called pythag:
function Complex_Abs (z: Complex): Double; var absx, absy: double; begin absx:=abs(Re(z)); absy:=abs(Im(z)); if absx>absy then Complex_Abs:=absx*sqrt(1.0+sqr(absy/absx)) else if absy=0 then Complex_Abs:=0.0 else Complex_Abs:=absy*sqrt(1.0+sqr(absx/absy)); end;
But then I've got to reverse Emil's optimization in Complex_Ln, so it won't fail for the same problem. Or there a better way to do complex Ln?
Now I had a look for this problem, both in the djgpp libc and libm libraries, where there is a function hyphot which does the job, and is said to be ieee754 standard, and in the CEPHES long double math library, (which defines cabs), which I use (with an import unit) when I need accuracy, since djgpp has no long double library (only elementary + - / * operations I think). libc overcomes the problem by using the extended accuracy and exponent range of the x87 floating point unit. gpc cannot use directly the same trick, since sqrt is computed in double accuracy only because djgpp does not define sqrtl (this is tested by configure). One way out would be to use an other configure test for hypot, to use directly this libc function (if this is possible).
Yes, I just hope it's implemented in a "good" way, otherwise it would be worse to use this function compared to the code above.
I'm doing both, and it works on my system. (maur8.pas)
Frank
-
CBFalconer -
Emil Jerabek -
Frank Heckenbach -
Maurice Lombardi