OK, if you're happy that the code is correct ... I guess I just don't understand what the "a := Abs(Im(z)) / (2 * b)" statement does ... More accurate results are always welcome, anyway. Anyway, although you might say the Clipped_SqRt solution is a kludge, it is certainly not second-guessing. It would only be applied where the minimum theoretical argument was zero, in which case it would avoid (and eliminate) rounding errors which produced a slightly negative "zero". Joe.

-----Original Message----- From: Frank Heckenbach [SMTP:frank@g-n-u.de] Sent: Thursday, January 10, 2002 9:30 PM To: gpc@gnu.de Subject: Re: Complex_SqRt bug

da Silva, Joe wrote:

...

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 ...

I don't think so. This seems more like a kludge to me (first compute a wrong value, and then second-guess in a subroutine what it could really be). Besides, it would only avoid the problem in the "< 0" case -- as Emil noted, "the patch above gives more accurate results (even for arguments which do not trigger the bug)."

BTW, I checked his code carefully to make sure it's mathematically correct.

Frank

-- Frank Heckenbach, frank@g-n-u.de, http://fjf.gnu.de/, 7977168E GPC To-Do list, latest features, fixed bugs: http://agnes.dida.physik.uni-essen.de/~gnu-pascal/todo.html