From owner-freebsd-numerics@FreeBSD.ORG Sun Aug 12 23:11:08 2012 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 7963D1065675 for ; Sun, 12 Aug 2012 23:11:08 +0000 (UTC) (envelope-from peter@rulingia.com) Received: from vps.rulingia.com (host-122-100-2-194.octopus.com.au [122.100.2.194]) by mx1.freebsd.org (Postfix) with ESMTP id F283C8FC15 for ; Sun, 12 Aug 2012 23:11:07 +0000 (UTC) Received: from server.rulingia.com (c220-239-249-137.belrs5.nsw.optusnet.com.au [220.239.249.137]) by vps.rulingia.com (8.14.5/8.14.5) with ESMTP id q7CNB7AP075866 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Mon, 13 Aug 2012 09:11:08 +1000 (EST) (envelope-from peter@rulingia.com) X-Bogosity: Ham, spamicity=0.000000 Received: from server.rulingia.com (localhost.rulingia.com [127.0.0.1]) by server.rulingia.com (8.14.5/8.14.5) with ESMTP id q7CNB1ZZ021700 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Mon, 13 Aug 2012 09:11:01 +1000 (EST) (envelope-from peter@server.rulingia.com) Received: (from peter@localhost) by server.rulingia.com (8.14.5/8.14.5/Submit) id q7CNB1Nl021699 for freebsd-numerics@freebsd.org; Mon, 13 Aug 2012 09:11:01 +1000 (EST) (envelope-from peter) Resent-From: Peter Jeremy Resent-Date: Mon, 13 Aug 2012 09:11:01 +1000 Resent-Message-ID: <20120812231101.GF20453@server.rulingia.com> Resent-To: freebsd-numerics@freebsd.org Received: from vps.rulingia.com (host-122-100-2-194.octopus.com.au [122.100.2.194]) by server.rulingia.com (8.14.5/8.14.5) with ESMTP id q6KIr4Oq058321 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Sat, 21 Jul 2012 04:53:04 +1000 (EST) (envelope-from brde@optusnet.com.au) Received: from mail26.syd.optusnet.com.au (mail26.syd.optusnet.com.au [211.29.133.167]) by vps.rulingia.com (8.14.5/8.14.5) with ESMTP id q6KIr49H094081 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Sat, 21 Jul 2012 04:53:04 +1000 (EST) (envelope-from brde@optusnet.com.au) 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 mail26.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q6KIqhvu007254 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Sat, 21 Jul 2012 04:52:45 +1000 From: Bruce Evans Mail-Followup-To: freebsd-numerics@freebsd.org X-X-Sender: bde@besplex.bde.org To: Stephen Montgomery-Smith In-Reply-To: <50097128.6030405@missouri.edu> Message-ID: <20120721032448.X5744@besplex.bde.org> References: <20120714120432.GA70706@server.rulingia.com> <20120717084457.U3890@besplex.bde.org> <5004A5C7.1040405@missouri.edu> <5004DEA9.1050001@missouri.edu> <20120717200931.U6624@besplex.bde.org> <5006D13D.2080702@missouri.edu> <20120718205625.GA409@troutmask.apl.washington.edu> <500725F2.7060603@missouri.edu> <20120719025345.GA1376@troutmask.apl.washington.edu> <50077987.1080307@missouri.edu> <20120719032706.GA1558@troutmask.apl.washington.edu> <5007826D.7060806@missouri.edu> <5007AD41.9070000@missouri.edu> <20120719205347.T2601@besplex.bde.org> <50084322.7020401@missouri.edu> <20120720035001.W4053@besplex.bde.org> <50085441.4090305@missouri.edu> <20120720162953.N2162@besplex.bde.org> <20120720184114.B2790@besplex.bde.org> <50097128.6030405@missouri.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Mailman-Approved-At: Sun, 12 Aug 2012 23:56:00 +0000 Cc: Diane Bruce , Bruce Evans , John Baldwin , David Chisnall , Bruce Evans , Steve Kargl , David Schultz , Peter Jeremy , Warner Losh Subject: Re: Use of C99 extra long double math functions after r236148 X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Date: Sun, 12 Aug 2012 23:11:08 -0000 X-Original-Date: Sat, 21 Jul 2012 04:52:43 +1000 (EST) X-List-Received-Date: Sun, 12 Aug 2012 23:11:08 -0000 On Fri, 20 Jul 2012, Stephen Montgomery-Smith wrote: > On 07/20/2012 04:19 AM, Bruce Evans wrote: >> % x0 = ldexp(floor(ldexp(x, 24)), -24); >> % x1 = x - x0; >> % y0 = ldexp(floor(ldexp(y, 24)), -24); >> % y1 = y - y0; >> >> This has a chance of working iff the bound away from 1 is something like >> 2**-24. Otherwise, multiplying by 2**24 and flooring a positive value >> will just produce 0. 2**-24 seems much too small a bound. My test >> coverage is not wide enough to hit many bad cases. > > This is meant to cover a situation where x = cos(t) and y = sin(t) for some t > not a multiple of PI/2. > > Now, hypot(x,y) will be 1, but only to within machine precision, i.e. an > error of about 1e-17. Actually more like DBL_MIN = 2e-308 for the real part (2e-321 for the smallest denormal). If you are sloppy and have an error of 1e-17, then the relative error is a 2e304 times :-). In ulps I only saw errors below 1e20 ulps. > So log(hypot(x,y)) will be about 1e-17. The true answer being 0, the ULP > will be infinite. The scaling for a correctly rounded result of 0 is unclear. But when the correctly rounded result is the smallest denormal, the scaling is clear: a result of 0 is off by 1 ulp, and a result of 1e-17 is off by 2e304 ulps :-). Both with a factor of 2 of fuzziness for the limited precision of the result. > BUT (and this goes with Goldberg's paper when he considers the quadratic > formula when the quadratic equation has nearly equal roots), suppose > > x = (double)(cos(t)) > > that is, x is not exactly cos(t), but it is the number "cos(t) written in > IEEE double precision". Similarly for y. That is, even though the formula > that produce x and y isn't exact, let's pretend that x and y are exact. > > Again log(hypot(x,y)) will be about 1e-17. But the true answer will also be > about 1e-17. But they won't be the same, and the ULP will be about 1e17. > > What my formula does is deal with the second case, but reduce the ULP to > about 1e8! That is, if x and y are exact numbers, and it so happens that > hypot(x,y) is very close to 1, my method will get you about 8 extra digits of > accuracy. It was still giving errors of 1e17 (maybe 1e304) ulps for me. > Now you have special formulas that handle the cases when z is close to 1, -1, > I and -I. That was my old version. Now I use your formula in most cases. It seems to give 300 or so extra digits after debugging it. I still need the special formula the smallest denormal result. That may be the only case (when y*y/2 == 0 instead of the smallest denormal). > Earlier this morning I sent you a formula, which I think might be > slightly more accurate than yours, for when z is close to 1. I think similar > formulas can be produced for when z is close to -1, I and -I. I think this only gives a small error relative to |log(z)|, so it gives a huge error relative to real(log(z)) when that is nearly 0. Indeed, the error increased from 1 ulp to 2e16 ulps. It can't handle the case of the smallest denormal: z = 1 + I*y where y*y/2 ~= smallest denormal. Then |z - 1| is |y| which is about 1e-162, and log(|z|) is about y*y/2 = smallest denormal. Denormals only go a few powers of 10 lower here. This magic 2e16 is essentially 2**(DBL_MANT_DIG + 1). Errors in ulps are never of the order of 1e304 or 1e162. > To get ULP of about 1 when x and y are exact, and it happens that hypot(x,y) > is close to 1, but z is not close to 1, -1, I or -I, would require, I think, > hypot(x,y)-1 being computed using double double precision (i.e. a mantissa of > 108 bits), and then feeding this into log1p. It does need about that, but this is routine for hi-lo decompositions, and you claimed to already have it :-) : @ x0 and y0 are good approximations to x and y, but have their bits trimmed @ so that double precision floating point is capable of calculating @ x0*x0 + y0*y0 - 1 exactly. */ Bruce