math function speed

[ark at alice.UUCP]

In article <2096 at ulysses.homer.nj.att.com>, ggs at ulysses.UUCP writes:
-> In article <5640 at brl-smoke.ARPA>, gwyn at brl-smoke.UUCP writes:
-> > 
-> > In a test where the result of "sqrt" was compared to known correct
-> > answers, for the following arguments:
-> > 
-> > 	0	1/5	1/3	1/2	1	2	e
-> > 	3	pi	4	5	9	10	16
-> > 	49	100	10000
-> > 
-> > I got the following results:
-> > 
-> > SVR2 function "sqrt":
-> > 
-> >     max abs error:	  2.7755576e-17
-> >     max rel error:	  2.4037034e-17
-> > 
-> > 4.3BSD function "sqrt":
-> > 
-> >     max abs error:	  5.5511151e-17
-> >     max rel error:	  3.3669215e-17
-> 
-> I thought the usual measure of accuracy was in units of "least
-> significant bits".  I repeated my test that compared a "best" solution
-> with a library version, this time with the sqrt that I got from the BRL
-> System V emulation package (sorry, I don't have a real VAX System V to
-> try, and the 3B20 is a different beast).  My result, using 100000 random
-> positive doubles:
-> 
-> SVR2 function "sqrt"
-> 
->     one bit errors:	20356
->     two bit errors:	0
-> 
-> 4.3BSD function "sqrt"
-> 
->     one bit errors:	8634
->     two bit errors;	0
-> 
-> I question your choice of test numbers: any sqrt function that misses
-> 0 is broken.  Seven of the seventeen numbers are perfect squares,
-> and six are powers of two.  They are excellent boundary value tests, but
-> statistically unusual.

The IEEE P754 definition of square root requires that the result be
the most accurate possible.  That is, what you get must be exactly
equal to the infinite-precision square root, correctly rounded.

This is not as difficult as it may sound, because the rounding can
never result in a tie.  That is, the infinite-precision square root
of a floating-point number can never fall exactly between two
adjacent floating-point numbers, so you only need one extra bit
of precision to get the right answer.

To see that this is true, consider Y = sqrt(X), after rounding, and
scale Y by a a power of 2 (and X by that power of 2 squared)
until the low-order bit of Y represents 1.  Now, assume that
the rounding of Y came after a tie.  In other words, in infinite
precision, X = (Y + 0.5) ** 2.  Can this be possible if both X
and Y are integers?  (that's why I did the scaling)

No, it can't be possible.  (Y + 0.5) ** 2 is Y**2 + 2*Y*0.5 + 0.25,
which can't be an integer.

So for floating-point square root, it is entirely meaningful
to ask whether or not the results are correct, and looking at
relative error just conceals the real issue.