From owner-freebsd-numerics@FreeBSD.ORG Sun Aug 12 23:11:04 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 1DDBC106564A for ; Sun, 12 Aug 2012 23:11:04 +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 7D37F8FC14 for ; Sun, 12 Aug 2012 23:11:03 +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 q7CNB3nK075850 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Mon, 13 Aug 2012 09:11:03 +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 q7CNAvAo021684 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Mon, 13 Aug 2012 09:10:57 +1000 (EST) (envelope-from peter@server.rulingia.com) Received: (from peter@localhost) by server.rulingia.com (8.14.5/8.14.5/Submit) id q7CNAv6H021683 for freebsd-numerics@freebsd.org; Mon, 13 Aug 2012 09:10:57 +1000 (EST) (envelope-from peter) Resent-From: Peter Jeremy Resent-Date: Mon, 13 Aug 2012 09:10:56 +1000 Resent-Message-ID: <20120812231056.GE20453@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 q6KEssvS031624 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Sat, 21 Jul 2012 00:54:55 +1000 (EST) (envelope-from stephen@missouri.edu) Received: from wilberforce.math.missouri.edu (wilberforce.math.missouri.edu [128.206.184.213]) by vps.rulingia.com (8.14.5/8.14.5) with ESMTP id q6KEsqnW086889 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Sat, 21 Jul 2012 00:54:54 +1000 (EST) (envelope-from stephen@missouri.edu) Received: from [127.0.0.1] (wilberforce.math.missouri.edu [128.206.184.213]) by wilberforce.math.missouri.edu (8.14.5/8.14.5) with ESMTP id q6KEsWie037880; Fri, 20 Jul 2012 09:54:32 -0500 (CDT) (envelope-from stephen@missouri.edu) Message-ID: <50097128.6030405@missouri.edu> From: Stephen Montgomery-Smith Mail-Followup-To: freebsd-numerics@freebsd.org User-Agent: Mozilla/5.0 (X11; Linux i686; rv:14.0) Gecko/20120714 Thunderbird/14.0 MIME-Version: 1.0 To: Bruce Evans 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> In-Reply-To: <20120720184114.B2790@besplex.bde.org> Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit Cc: Diane Bruce , 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:04 -0000 X-Original-Date: Fri, 20 Jul 2012 09:54:32 -0500 X-List-Received-Date: Sun, 12 Aug 2012 23:11:04 -0000 On 07/20/2012 04:19 AM, Bruce Evans wrote: > On Fri, 20 Jul 2012, Bruce Evans wrote: > The only way for x*x + y*y to be _very_ near 1 in infinite precision > is for |x| or y| to be 1 (I think). Other cases are bounded away from > 1, and if you are lucky the bound is fairly far from 1 so that sloppier > approximations work OK. Mathematicians should determine the bound > exactly using continued fractions or something like they do for > approximations to N*Pi/2. This becomes especially interesting in high > precisions where you can't hope to get near the worst case by random > testing. > > % 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. So log(hypot(x,y)) will be about 1e-17. The true answer being 0, the ULP will be infinite. 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. Now you have special formulas that handle the cases when z is close to 1, -1, I and -I. 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. 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. Of course, you may argue that situations when x and y are exact, not close to 1 or 0, and hypot(x,y) is close to 1, are so very rare that extra consideration is not required. But my algorithm produces better answers than the naive formula even when the distance between 1 and hypot is about 1/10. The naive formula has a ULP of about 10, and I get it down to less than 2. And when the distance between hypot(x,y) and 1 is about 1e-5, the naive formula has a ULP of about 1e5, and I still manage to get a ULP of about less than 2. Stephen