Re: libm-test.inc: Computing ulps near FP_ZERO.

["Carlos O'Donell" <carlos at redhat dot com>]

On 04/06/2013 06:42 PM, Brooks Moses wrote:
> Carlos O'Donell wrote, at 4/6/2013 11:46 AM:
>> testing double (inline functions)
>> Failure: Test: cos (pi/2) == 0
>> Result:
>>  is:          6.12323399573676603587e-17   0x1.1a62633145c070000000p-54
>>  should be:   0.00000000000000000000e+00   0x0.00000000000000000000p+0
>>  difference:  6.12323399573676603587e-17   0x1.1a62633145c070000000p-54
> [...]
>> As expected 6.123e-17 is a *terrible* answer for zero. It should be
>> a failure until we get down to < 16 ulp.
> [...]
>> Did I miss something or misunderstand some concept?
> 
> You've missed something.  :)
> 
> What you've missed here is that the "pi/2" input is not equal to the
> mathematical pi/2, because the latter is not representable as a
> floating-point number.  The deviation is bounded by ulp(pi/2)/2, and
> since the first derivative of cosine at pi/2 is 1, that means that an
> _exact_ cosine function will produce a result that is off of 0 by an
> amount that is likewise bounded by ulp(pi/2)/2.

This is a very good point, one which I wasn't thinking about as
I worked through the problem of computing ulps in glibc.

In this particular case I would argue the test suite is wrong and
instead of expecting zero it should expect ulp(pi/2)/2, given
that one of the expected input arguments is impossible to represent
exactly. That makes sense.

It still doesn't answer the question that I'm trying to raise here,
which is how do we compute ulps near FP_ZERO?

Take another example, particularly the one from the LSB oliver
testsuite which does w^z, where w = -1 -0*i and z = -1 -0*i.
In this case both inputs are perfectly representable as floating
point arguments but the answer is significantly off from the
expected zero. In this case I can make no easy assumptions that
the input arguments dictate that I adjust the answer away from
-1 + 0*i in any appreciable way.

The pen-and-paper answer to this is as follows:

w = -1 -0*i
t = arg(w) = pi
r = |w| = 1
z = c + d*i
z = -1 -0*i
c = -1
d = -0
w^z = e^(z*log(w))
w^z = e^(z*(log(r)+t*i)
w^z = r^(c)*e^(-d*i)[cos(d*log(r)+c*t)+i*sin(d*log(r)+c*t)]
w^z = 1^(-1)*e^(-0*i)[cos(-0*log(1)-1*pi)+i*sin(-0*log(1)-1*pi)]
w^z = 1*1*[cos(-0-pi)+i*sin(-0-pi)]
w^z = -1 + 0*i

Were I to have used this solution to compute w^z one can immediately
see that we face the same problem as cos(pi/2), in that pi is not
representable and thus the result of sin(-0-pi) is ulp(pi)/2
away from zero or ~2.2e-16 (0x1.0p-52).

I can't know this apriori though. The algorithm might change,
and several other factors come into play here.

In the case of cos(pi/2) the answer seems clear: the test case must
take into account that the input argument is not precisely
representable.

In the case of w^z, I see no way to make that argument. I must
request that the test case expect -1 + 0*i, and in that case I must
be able to compare ~0*i to +0*i and give some ulp value for this.

That is to say: How many ulps away from +0*i is ~0*i?

We are back to this:

java says ulp(0x0.0p0) is 0x0.0p-1022, which is the next step to
a normal fp value.

I expected ulp(0x0.0p0) to be 0x0.0p-1074, which is the next step
to the smallest subnormal fp value.

glibc says ulp(0x0.0p0) is 0x1.0p-52, which is a huge ulp value to
use around zero.

Is glibc's choice of ulp(0x0.0p0) in order to account for the fact
that the test suite expects 0x0.0p0 from cos (pi/2), when that 
isn't true (or even possible)?

Cheers,
Carlos.