Date: Sun, 29 Jul 2012 06:10:03 GMT From: Bruce Evans <brde@optusnet.com.au> To: freebsd-bugs@FreeBSD.org Subject: Re: bin/170206: complex arcsinh, log, etc. Message-ID: <201207290610.q6T6A3lV085252@freefall.freebsd.org>
next in thread | raw e-mail | index | archive | help
The following reply was made to PR bin/170206; it has been noted by GNATS.
From: Bruce Evans <brde@optusnet.com.au>
To: Stephen Montgomery-Smith <stephen@missouri.edu>
Cc: Bruce Evans <brde@optusnet.com.au>, freebsd-bugs@FreeBSD.org,
FreeBSD-gnats-submit@FreeBSD.org,
Stephen Montgomery-Smith <stephen@FreeBSD.org>
Subject: Re: bin/170206: complex arcsinh, log, etc.
Date: Sun, 29 Jul 2012 15:59:49 +1000 (EST)
On Sat, 28 Jul 2012, Stephen Montgomery-Smith wrote:
> On 07/28/2012 11:15 AM, Stephen Montgomery-Smith wrote:
>> On 07/28/2012 10:46 AM, Stephen Montgomery-Smith wrote:
>>> OK. This clog really seems to work.
>>>
>>> x*x + y*y - 1 is computed with a ULP less than 0.8. The rest of the
>>> errors seem to be due to the implementation of log1p. The ULP of the
>>> final answer seems to be never bigger than a little over 2.
> ...
> And I should add that I just realized that ULP isn't quite the same as
> relative error, so an extra factor of up to 2 could make its way into the
> calculations.
Yes, this is tricky. For denormals, it is easy to be off by a factor of
2**MANT_DIG or infinity instead of only 2.
For normals, the most interesting cases are near powers of 2 (say 1).
One ulp is twice as large for values in [1, 2) as it is for values
in [0.5, 1). Even to determine which one to use, you need to know
if the infinitely precise result is >= 1 or < 1, else you may be
off by a factor of 2 in the error checking. If the factor is 1/2,
then it hides errors, and if it is 2 then it gives unexpected errors.
For denormals, the easiest case to understand is when the correctly
rounded case is the smallest strictly positive denormal. Then the
size of an ulp is the same as the value of this denormal. A rounding
error of < 1 ulp (but > 0.5 ulps) may give a result of 0 ulps or 2 ulps.
Such errors are to be expected. But relatively, they are infinity and
2, respectively. Normally you expected rounding errors of near
2**-MANT_DIG, and infinity and 2 are much larger. The relative error
should be scaled (like the size of an ulp should be) so that you don't
see such large errors.
You might not care about denormals, but you should check a few of them
and then you don't want errors in the checking software for them hiding
errors for normals. Without denormals, there would be no gradual
underflow, and underflow from the smallest normal to 0 really would be
essentially infinity (it would be about 2**MANT_DIG in ulps).
My checking software scales for denormals, but I think I got it wrong
and am off by a factor of 2. For normals near a power of 2, its just
impossible to determine the right scale without a reference function
with enough extra precision to determine which side it is on. Even
the correct definition of an ulp is unclear near powers of 2 (perhaps
other cases). I want my checking program to match my definition, which
is currently just what the checking program does :-) -- don't worry
about the reference function not be precise enough. There is the
same uncertainty about the size of an ulp for a result of 0 (and
more -- maybe the result should be -0 :-).
Recently I started checking extra internal precision for some functions.
This gives relative errors of < 2**(-MANT_DIG+N), where N is the extra
precision. When the relative error is near 2**-MANT_DIG, it is hard
to tell if it is smaller than 1 ulp, since the ulp scale may be hard
to determine and the relative error doesn't map simply to ulps. When
N >= 2, there is a factor of 2 to spare and final errors (after
rounding) that are certain to be < 1 ulp are easy to ignore. Even
functions that don;t have much extra internal precision but which
that meet my design goal of a final error < 1 ulp should have N at
least 1, so their non-errors should be easy to not see too, but I have
used this mainly for functions with N about 7. The internal errors
for the recently committed logl() are < 2**-120 relative on sparc64,
except for denormal results. This is obviously enough for 113-bit
precision to 0.5+epsilon ulps. But rounding would cluster the final
errors near 2**-113 or 2**--114 and it wouldn't be obvious how this
maps to ulps.
Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?201207290610.q6T6A3lV085252>