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