math function speed
ark at alice.UUCP
ark at alice.UUCP
Sun Feb 22 02:35:42 AEST 1987
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.
More information about the Comp.lang.c
mailing list