Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Path: utzoo!watmath!clyde!rutgers!princeton!allegra!alice!ark
From: ark@alice.UUCP
Newsgroups: comp.lang.c
Subject: Re: math function speed
Message-ID: <6655@alice.uUCp>
Date: Sat, 21-Feb-87 11:35:42 EST
Article-I.D.: alice.6655
Posted: Sat Feb 21 11:35:42 1987
Date-Received: Sun, 22-Feb-87 02:04:37 EST
References: <4943@mimsy.UUCP> <2550005@hpisod2.HP> <756@unc.unc.UUCP> <2096@ulysses.homer.nj.att.com>
Distribution: na
Organization: AT&T Bell Laboratories, Liberty Corner NJ
Lines: 67

In article <2096@ulysses.homer.nj.att.com>, ggs@ulysses.UUCP writes:
-> In article <5640@brl-smoke.ARPA>, gwyn@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.