Date: Sun, 12 Aug 2012 23:11:08 -0000 From: Bruce Evans <brde@optusnet.com.au> To: Stephen Montgomery-Smith <stephen@missouri.edu> Cc: Diane Bruce <db@db.net>, Bruce Evans <brde@optusnet.com.au>, 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: <20120721032448.X5744@besplex.bde.org> Resent-Message-ID: <20120812231101.GF20453@server.rulingia.com> In-Reply-To: <50097128.6030405@missouri.edu> 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>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120721032448.X5744>