This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Re: libm-test.inc: Computing ulps near FP_ZERO.
- From: Brooks Moses <brooks at codesourcery dot com>
- To: Carlos O'Donell <carlos at redhat dot com>
- Cc: GNU C Library <libc-alpha at sourceware dot org>, "Joseph S. Myers" <joseph at codesourcery dot com>, Andreas Jaeger <aj at suse dot com>, Thomas Schwinge <thomas at codesourcery dot com>, David Miller <davem at davemloft dot net>
- Date: Sat, 6 Apr 2013 15:42:31 -0700
- Subject: Re: libm-test.inc: Computing ulps near FP_ZERO.
- References: <51606D8A dot 9080003 at redhat dot com>
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.
(Since glibc's cosine function is representing its answer as a
floating-point number, we expect an additional error on the output equal
to ulp(output), which is negligible by comparison.)
Coincidentally, ulp(pi/2)/2 is exactly equal to ulp(1), which is what
you note that glibc is using in this case.
- Brooks