Date: Sun, 12 Aug 2012 23:11:04 -0000 From: Stephen Montgomery-Smith <stephen@missouri.edu> To: Bruce Evans <brde@optusnet.com.au> Cc: Diane Bruce <db@db.net>, John Baldwin <jhb@freebsd.org>, David Chisnall <theraven@freebsd.org>, Bruce Evans <bde@freebsd.org>, Steve Kargl <sgk@troutmask.apl.washington.edu>, David Schultz <das@freebsd.org>, Peter Jeremy <peter@rulingia.com>, Warner Losh <imp@bsdimp.com> Subject: Re: Use of C99 extra long double math functions after r236148 Message-ID: <50097128.6030405@missouri.edu> Resent-Message-ID: <20120812231056.GE20453@server.rulingia.com> In-Reply-To: <20120720184114.B2790@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>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?50097128.6030405>