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

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

On 04/08/2013 10:09 PM, Rich Felker wrote:
> 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.

Agreed. I was trying to avoid needing to use external software to
validate the result and trying to achieve this from a first principles
perspective.

The following patch fixes the cos failure.
---
@@ -5330,7 +5365,7 @@ cos_test (void)
 
   TEST_f_f (cos, M_PI_6l * 2.0, 0.5);
   TEST_f_f (cos, M_PI_6l * 4.0, -0.5);
-  TEST_f_f (cos, M_PI_2l, 0);
+  TEST_f_f (cos, M_PI_2l, 6.123233995736765886130329661375001e-17L);
 
   TEST_f_f (cos, 0.75L, 0.731688868873820886311838753000084544L);
 
---

The point that Brooks was making was that cos(pi/2) has an error,
namely that pi/2 is not equal to the floating point value we use to
represent pi/2, which means the output is never zero. The first
derivative of cos around pi/2 is ~1 and thus I was trying to see if
I could say the answer was 0+ulp(pi/2)/2.

To go back to first principles:

M_PIl: 
	3.1415926535897932384626433832795029L
M_PI_2l:
	3.1415926535897932384626433832795029L/2
	1.57079632679489661923132169163975145
double M_PI_2l:	
	1.5707963267948965579989817342720925807952880859375
	0x1.921fb54442d18p+1

The difference is:
          1.57079632679489661923132169163975145...
        - 1.5707963267948965579989817342720925807952880859375000...
	= 0.00000000000000006123233995736760320

Thus the answer should be 6.123233995736760320e-17 (almost your
answer with Wolfram Alpha).	

So cos (0x1.921fb54442d18p+1) ~= 6.123233995736760320e-17.

The error in pi/2 is ~ulp(pi/2)/2 or 1.11022302462515654042e-16.

The problem is then in using the pi/2 error as the answer is impossible
since it is only an upper bound on the error.

> 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.

The test as it is today pretends there is an actual "pi" and that's
part of the problem.

Cheers,
Carlos.