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

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

On 04/07/2013 08:28 PM, Brooks Moses wrote:
> Carlos O'Donell wrote, at 4/7/2013 3:58 PM:
>> Excellent discussion. 
> 
> Thanks!  Glad it was helpful.
> 
>> At this point I'm in agreement with you.
>>
>> The likely course of action is to do the following:
>>
>> (a) Adjust tests that rely on accurate input to generate accurate
>>     output, such as the cos(pi/2) test. These tests will have their
>>     expected result adjusted by the error involved in converting
>>     their input argument.
>>
>> (b) Enhance algorithms to detect zero branches. In particular fix
>>     cpow to detect zero branches. We have a microbenchmark framework
>>     now so I could run the new algorithms through that to see if
>>     the function is faster or slower.
>>
>> Comments?
> 
> FWIW, that sounds reasonable to me, but I'm just the peanut gallery.  :)
> 
> I'd be interested in knowing how the microbenchmarks come out.  My guess
> is that adding check for zero-branches (which you'd implement just by
> checking for imag(z)==0, and then do other stuff only if that's true)
> will be negligible cost on most architectures with branch prediction.
> And my suggested alternate implementation might be faster, because
> you're avoiding some trigonometric calls.

I'll post results as further patches.

>> What about the cases where we can't easily fix the algorithm to
>> detect the zero branches? What do we use for ulp(0x0.0p0) then? :-)
> 
> I figure those fall into three categories:
> 
>  * You're getting zero through something that amounts to subtracting
>    two "equal" quantities (in the above cases, that was a pi/2 input
>    and a pi/2 used to condition the input, basically), and the
>    expected error is based on ulp(input), not ulp(output).
> 
>  * A zero input propagates through to a zero output.  Here, you ought
>    to get exactly zero in the output.  This is what happened in the
>    algorithm I wrote out for w^z; it doesn't actually detect zeros in
>    a branching sense, but will preserve them nonetheless.
> 
>  * A zero output results from underflow.  Here you also ought to get
>    exactly zero.
> 
> Are there any cases in the testsuite that are exceptions to this
> categorization?

I don't think so.

> I think anything where you can't get an exact zero is going to end up
> being like the w^z-with-the-current-implementation case -- where the
> error is inherently on the order of ulp(1) or so, and the best the test
> can do is document why -- rather than giving you an error on the order
> of the smallest floating-point numbers.

The following is the only testsuite failures I see after fixing up
the definition of ulp to match what is expected from a mathematical
perspective.

It reveals quite clearly, that cos, sincos, and cpow have some
terrible errors.

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
Maximal error of `cos'
 is      : 2019720827359741 ulp
 accepted: 2 ulp
Failure: Test: sincos (pi/2, &sin_res, &cos_res) puts PI_ERROR/2 in cos_res
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
Maximal error of `sincos'
 is      : 2019720827359741 ulp
 accepted: 1 ulp
Failure: Test: Imaginary part of: cpow (e + 0 i, 0 + 2 * M_PIl i) == 1.0 - PI_ERROR i
Result:
 is:         -2.44929359829470641435e-16  -0x1.1a62633145c070000000p-52
 should be:  -2.22044604925031308085e-16  -0x1.00000000000000000000p-52
 difference:  2.28847549044393333500e-17   0x1.a62633145c0700000000p-56
 ulp(x)    :  4.93038065763132378382e-32   0x1.00000000000000000000p-104
 ulp       :  464157972651015.0000
 max.ulp   :  0.0000
Maximal error of real part of: cpow
 is      : 2 ulp
 accepted: 2 ulp
Maximal error of imaginary part of: cpow
 is      : 464157972651015 ulp
 accepted: 2 ulp

Test suite completed:
  7881 test cases plus 7253 tests for exception flags executed.
  6 errors occurred.


Cheers,
Carlos.

Attachment: libm-test.inc.diff-v4
Description: Text document