Skip site navigation (1)Skip section navigation (2)
Date:      Sun, 29 Jul 2012 02:20:03 GMT
From:      Bruce Evans <brde@optusnet.com.au>
To:        freebsd-bugs@FreeBSD.org
Subject:   Re: bin/170206: complex arcsinh, log, etc.
Message-ID:  <201207290220.q6T2K3an059742@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 12:14:18 +1000 (EST)

 On Sat, 28 Jul 2012, 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.
 >
 > Also, I don't think the problem is due to the implementation of log1p. If you 
 > do an error analysis of log(1+x) where x is about exp(-1)-1, and x is correct 
 > to within 0.8 ULP, I suspect that about 2.5 ULP is the best you can do for 
 > the final answer:
 >
 > relative_error(log(1+x)) = fabs(1/((1+x) log(1+x))) * relative_error(x)
 >                  = 1.58 * relative_error(x)
 >
 > Given that log1p has itself a ULP of about 1, and relative error in x is 0.8, 
 > and considering x=exp(-1)-1, this gives a ULP at around 1.58*0.8+1 = 2.3. 
 > And that is what I observed.
 
 Not given:
 - the error in log1p is close to 0.5 ulps in my version
 - when the extra internal precision of my log1p is brought out, the error
    will be a maximum of about 1/2**7 ulps before rounding.
 This will only save about half an ulp after rounding.  The are accumalative
 error of about half an ulp from the other 2+ steps).  Extra precision for
 these is further off.
 
 I made progress towards completing exhaustive testing of this point
 for all cases near 1 for clogf.  All cases tested (at least 99% of
 all cases possible) were perfectly rounded.  There were an amazingly
 (before doing the analysis below) large number of half-way cases
 near (1 + ~1e-10*I).  This is an easy special case for clogf() -- x
 is precisely 1 so it reduces immediately to log1pf(1+y*y)/2 -- but
 when x is precisely 1 there are many half-way cases.  Due to
 numerical accidents, it turns out that many cases are correctly
 rounded in float precision but not in double precision.  My testing
 is not quite complete because it doesn't verify that the accidents
 aren't mostly due to my test methodology.  Part of the methodology
 is:
 
 - start with float args.  These have only 24 mantissa bits, so after
    promotion to double and long double precision they have many lower bits
    0.  This is what gives they many half-way cases that are unresolvable
    in long double precision after y is squared.  The magic arg values
    (y near 1e-10) are what is needed for squaring to push the lowest 1
    bit into an unfortunate place.
 - in both clogf() and clog(), return log1pl((long double)y*y)/2 (rounded
    to float precision), so that we see the inherent inaccuracy of clogf()
    and not external inaccuracy from log1pf().  This is uninteresting
    for x = 1 since there are no internal extra-precision calculations
    to verify for that case, but it is interesting for the general case
    to verify that the doubled float precision algorithm is exact.  I
    don't want the general case hidden by errors for this case, so I
    evaluate y*y in long double precision so as to minimise rounding
    errors in it, although the normal inaccurate y*y is part of clogf()
    (we could use doubled precision for it, but we don't because we know
    that all the extra precision would be lost when it is passed to
    log1pf().  Except when y*y is near underflow, we are more careful).
 - it turns out that the extra precision of log1pl() is enough for
    perfect rounding in float precision but not for perfect rounding
    in double precision.  This is no accident.  log1pl() gives only
    11 extra bits for double precision, but 40 extra for float
    precision.  A simple probabistic argument shows that 40 is always
    enough unless we are unlucky.  The chance of a half-way case that
    is not resolvable in long double precision is about 2*--11 and
    2**-40 in the 2 cases, respectively.  There are less than 2**-24
    args of interest (all x between sqrt(2)-epsilon and 1, corresponding
    y near sqrt(1-x*x)).  We would have to be much more unlucky to hit
    a half-way case with so many fewer cases in float precision.
    Assuming that the half-way cases are uniformly distributed, which
    they aren't, then the probabities for _not_ hitting a half-way case
    would be related to:
      (1-2**-40)**(2**24) = 0.999985  for float precision
      (1-2**-11)**(2**53) = <too small for pari>, but much smaller than
      (1-2**-11)**(2**39) = 3.315e-116608509  for double precision
    This calculation is certainly off by many powers of 2 even in float
    precision.  It's a bit surprising that the probability is so high
    for 2**24 cases (no birthday paradox with uniform distribution).
 
 Bruce



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?201207290220.q6T2K3an059742>