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
da Silva, Joe wrote:
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 ...
I could write down the computations if you like (unless Emil has them readily available already) ...
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".
As I said, I'm no expert on numeric algorithms, but I think if an algorithm can produce such wrong results, it might also be inaccurate for other values and should be improved (such as Complex_Sqrt). After all, normal rounding doesn't easily make a positive value (even very close to 0) negative. So it might be better if the programmer sees the problem clearly with an error message. If it really turns out that an algorithm is as good as possible (numerically) and still prone to this kind of problem, something like Clipped_SqRt might be in order. But I guess this is rather an exceptional situation, so Clipped_SqRt does not seem a good candidate for a general library routine.
But perhaps Emil can comment better on these matters ...
Frank
Joe da Silva wrote:
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 ...
That's easy. Let u+iv=\sqrt{x+iy}. Then (u+iv)^2=x+iy, i.e. y=2uv. (The code makes either a:=u and b:=v, or a:=v and b:=u, depending on the sign of x. Both "a := ..." assignments in the "if" statement are mathematically equivalent, but "SqRt ((a - Abs (Re (z))) / 2)" is too inexact if |y|<<|x|, since it computes a _small_ number as a difference of two _big_ numbers.)
Emil
On Mon, 14 Jan 2002, Emil Jerabek wrote:
Joe da Silva wrote:
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 ...
That's easy. Let u+iv=\sqrt{x+iy}. Then (u+iv)^2=x+iy, i.e. y=2uv. (The code makes either a:=u and b:=v, or a:=v and b:=u, depending on the sign of x. Both "a := ..." assignments in the "if" statement are mathematically equivalent, but "SqRt ((a - Abs (Re (z))) / 2)" is too inexact if |y|<<|x|, since it computes a _small_ number as a difference of two _big_ numbers.)
I guess this is because with 16-digit precision double floating-point numbers things like
(1e20 + 1 - 1 == 1e20)
fail ....
Also,
(1 + 0.1 + ... + 1/(10^n) == 1/(10^n) + ... + 0.1 + 1)
... can fail - the right hand side calculation is more precise if both are done left-to-right.
But it can be that GNU libm on Linux is doing it in the order shown on left hand side, since by using 'long double' intermediate results and adding Taylor series for cos(x) I have achieved _exactly the same_ square error sum of
sum (x - acos(cos(x )))^2 i=1..n i i
with my Taylor series as with original libm's cos(x)!!!!
If the sums of squares of errors are exactly the same over a wide range of parameters, my guess is that they're using the same Taylor series, and that they're adding it from larger to smaller series member -- yes, as you guessed: IN LESS PRECISE ORDER!
Now, could it be that GNU libm on x86-*-linux is just using Intel FPU's result, and that FPU on x86 chips is suffering from this imprecision?
Any ideas?
Mirsad
-- This message has been made up using recycled ideas and language constructs. No tree has been cut nor animal harmed in process of making it.
Visit http://www.hungersite.com/, save one CHILD from DYING!!!
-
da Silva, Joe -
Emil Jerabek -
Frank Heckenbach -
Mirsad Todorovac