bin/170206: complex arcsinh, log, etc.

[Bruce Evans]

The following reply was made to PR bin/170206; it has been noted by GNATS.

From: Bruce Evans <brde at optusnet.com.au>
To: Stephen Montgomery-Smith <stephen at missouri.edu>
Cc: Bruce Evans <brde at optusnet.com.au>, freebsd-bugs at FreeBSD.org,
        FreeBSD-gnats-submit at FreeBSD.org,
        Stephen Montgomery-Smith <stephen at 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