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

[Rich Felker <dalias at aerifal dot cx>]

On Mon, Apr 08, 2013 at 06:38:23PM -0700, Brooks Moses wrote:
> Brooks Moses wrote, at 4/8/2013 6:04 PM:
> > Carlos O'Donell wrote, at 4/8/2013 4:02 PM:
> >> Regenerating ULPs for /home/carlos/build/glibc/math/test-idouble
> >> testing double (inline functions)
> >> Failure: Test: cos (pi/2) == PI_ERROR/2
> >> Result:
> >>  is:          6.12323399573676603587e-17   0x1.1a62633145c070000000p-54
> >>  should be:   1.11022302462515654042e-16   0x1.00000000000000000000p-53
> >>  difference:  4.97899625051479936837e-17   0x1.cb3b399d747f20000000p-55
> >>  ulp(x)    :  2.46519032881566189191e-32   0x1.00000000000000000000p-105
> >>  ulp       :  2019720827359740.5000
> >>  max.ulp   :  0.0000
> > 
> > I think you've missed the bit where I pointed out that the expected
> > result here is "0 +/- ulp(pi/2)/2", not "ulp(pi/2)/2"!
> > 
> > In other words, you've confused the expected error bound with the
> > expected result.
> 
> Which is to say that there's a more fundamental problem here, because
> there's no clear way to implement what you actually need to do in the
> test harness as it stands.  This test suite is assuming (in, e.g.,
> check_float_internal) that some fixed multiple of ulp(expected) is a
> reasonable measure of the expected error bounds.
> 
> That, in this case, is fundamentally _false_.
> 
> You cannot meaningfully test cos(pi/2) with a function that takes the
> error and divides it by ulp(result) and then looks at the resulting
> number.  The result of that division is meaningless, because ulp(result)
> is irrelevant to the expected error.
> 
> Given those circumstances, setting ulp(0) equal to ulp(1) is a very
> workable hack around that limitation of the test harness, and will cause
> it to get the correct answer in the end -- namely, that a result of
> 6e-17 should be a pass, because ulp(1) is what we _actually_ want to be
> using as a threshold.
> 
> And, for the cases where the expected result is zero and you really do
> want exactly zero, you can just pass a zero max_ulp value to check_float
> and it will do what you need there too.

I'm a bit confused about this whole discussion. If the test is
cos(M_PI_2) then according to Wolfram Alpha,
cos(1.5707963267948965579989817342720925807952880859375), where that
argument is the exact value of M_PI_2, evaluates to approximately:

    6.123233995736765886130329661375001... Ã 10^-17

or, to the nearest double, 0x1.1a62633145c07p-54, which is greater
than the exact value. Thus, the result should be either
0x1.1a62633145c07p-54 or 0x1.1a62633145c06p-54, both of which are
within 1ulp of the correct value. Anything else is a buggy
implementation.

If on the other hand the test is pretending there's an actual "pi" to
work with, then the test itself is silly and should be replaced. It's
possible that you can reason about the error bounds for a conforming
implementation of cos based on the fact that the input is close to
pi/2 and using facts about the derivative, but this is overcomplicated
for a test.

Rich