From owner-freebsd-bugs@FreeBSD.ORG Sun Jul 29 02:14:28 2012 Return-Path: Delivered-To: freebsd-bugs@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [69.147.83.52]) by hub.freebsd.org (Postfix) with ESMTP id D1152106566C; Sun, 29 Jul 2012 02:14:28 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail04.syd.optusnet.com.au (mail04.syd.optusnet.com.au [211.29.132.185]) by mx1.freebsd.org (Postfix) with ESMTP id 618268FC08; Sun, 29 Jul 2012 02:14:28 +0000 (UTC) Received: from c122-106-171-246.carlnfd1.nsw.optusnet.com.au (c122-106-171-246.carlnfd1.nsw.optusnet.com.au [122.106.171.246]) by mail04.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q6T2EIVw012986 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Sun, 29 Jul 2012 12:14:20 +1000 Date: Sun, 29 Jul 2012 12:14:18 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Stephen Montgomery-Smith In-Reply-To: <50141018.3040203@missouri.edu> Message-ID: <20120729110356.Q1008@besplex.bde.org> 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> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-bugs@FreeBSD.org, FreeBSD-gnats-submit@FreeBSD.org, Stephen Montgomery-Smith , Bruce Evans Subject: Re: bin/170206: complex arcsinh, log, etc. X-BeenThere: freebsd-bugs@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Bug reports List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sun, 29 Jul 2012 02:14:28 -0000 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) = , 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