Skip site navigation (1)Skip section navigation (2)
Date:      Sun, 29 Jul 2012 12:14:18 +1000 (EST)
From:      Bruce Evans <brde@optusnet.com.au>
To:        Stephen Montgomery-Smith <stephen@missouri.edu>
Cc:        freebsd-bugs@FreeBSD.org, FreeBSD-gnats-submit@FreeBSD.org, Stephen Montgomery-Smith <stephen@FreeBSD.org>, Bruce Evans <brde@optusnet.com.au>
Subject:   Re: bin/170206: complex arcsinh, log, etc.
Message-ID:  <20120729110356.Q1008@besplex.bde.org>
In-Reply-To: <50141018.3040203@missouri.edu>
References:  <201207270247.q6R2lkeR021134@wilberforce.math.missouri.edu> <20120727233939.A7820@besplex.bde.org> <5012A96E.9090400@missouri.edu> <20120728142915.K909@besplex.bde.org> <50137C24.1060004@missouri.edu> <20120728171345.T1911@besplex.bde.org> <50140956.1030603@missouri.edu> <50141018.3040203@missouri.edu>

next in thread | previous in thread | raw e-mail | index | archive | help
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?20120729110356.Q1008>