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