Date: Fri, 7 Mar 2008 17:22:20 +1100 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Colin Percival <cperciva@FreeBSD.org> Cc: cvs-src@FreeBSD.org, src-committers@FreeBSD.org, Bruce Evans <bde@FreeBSD.org>, cvs-all@FreeBSD.org Subject: Re: cvs commit: src/sys/i386/include _types.h Message-ID: <20080307163246.J22709@besplex.bde.org> In-Reply-To: <47D089A4.4060208@freebsd.org> References: <200803051121.m25BLE03035426@repoman.freebsd.org> <20080305182531.GS68971@server.vk2pj.dyndns.org> <20080306021222.GA46783@zim.MIT.EDU> <47CF5D19.3090100@freebsd.org> <20080306033246.GA47280@zim.MIT.EDU> <47CF7EBF.6000009@freebsd.org> <20080306063452.GB48339@zim.MIT.EDU> <47CF9586.70707@freebsd.org> <20080307085031.P10724@delplex.bde.org> <47D089A4.4060208@freebsd.org>
next in thread | previous in thread | raw e-mail | index | archive | help
On Thu, 6 Mar 2008, Colin Percival wrote: > Bruce Evans wrote: >> On Wed, 5 Mar 2008, Colin Percival wrote: >>> You could have stopped this sentence here -- for all practical purposes, >>> correctly rounded trigonometric functions are not feasible. >> >> Nah, it's quite feasible, especially for transcendental functions since >> transcendental numbers in general can be approximated more closely by >> rationals than algebraic numbers. > > In fact, one of the common ways to prove that a real number is transcendental > is to show that it can be approximated more closely by rationals than any > algebraic function can be. I've only studied a few pages related to this. >> The evaluation just needs to be >> done in extra precision > > How much extra precision do you need? If the exact value of a transcendental > function is extremely close to halfway between two consecutive floating-point > values, how do you decide which way to round it? Probably not much more than the same number of extra bits as there are in the original precision. Values of trig functions mod 2^N are apparently sort of uniformly distributed. Thus each extra bit of precision gives correct rounding for about half of the remaining cases and N extra bits gives correct rounding for about 2^N cases, which is all cases if there were only 2^N cases to start with. This works out in practice for float precision -- with about 24 bits of extra precision, all except about 10 cases are perfectly rounded, and it is easy to evaluate those cases in extra precision to see that only a few more bits are needed. If the exact value is very close to half-way then you may need many more bits but would be unlucky to need many more. AFAIK (not far), no one knows exactly how many bits are needed for the worst case for even double precision. The domain is just too large to search (almost 2^64 bits). However, the large size of the domain means that very-near- half-way cases are very unlikely to occur in practice, so calling a slow multi-precision library to handle such cases would be a good way to handle them and sending mail about them would be a not so good way. But detecting such cases is too costly to do routinely, at least without hardware support. > Computing transcendental functions to within (0.5+x)ulp for small positive > values of x is certainly feasible, but computing provably correctly rounded > transcendental functions is very very hard. I think it is only very hard globally and for general functions. Elementary transcendental functions like exp() have rational coefficients and well-known properties which make them easier to bound. >>> Fortunately, library functions aren't required to have any particular >>> error bounds. >> >> We require error bounds of < 1 ulp and aim for < 0.5[0]1 ulps. > > What standard states that those bounds are required (or can be relied upon > when proving that code is correct)? Just the principle of least error. > Intel FPUs over the years have been > variously documented to compute transcendental functions to within 1.0, 2.5, > or 3.5 ulp -- are you saying that the 8087 doesn't comply with the standard > which it inspired? I didn't say that before, but in fact the ix87 has always had multi-gulp errors where its docs indicate errors of < 1 ulp, mainly for trig functions. Its trig functions can't possibly have errors of < 1 ulp, since its approximation to pi has only 66-68 bits, but an approximation with thousands of bits is needed (see libm). Multi-gulp errors occur as early as for x near pi/2 in 64-bit precision, since subtracting a 64-bit x near pi/2 from a 66 or 68 bit approx to pi/2 gives about 2 or 4 bits of accuracy or about 62 or 60 bits of error. 2^60 is about 1 exa-ulp (10^9 gulps). The errors are much larger for larger x, and give multi-gulp errors even for float precision for not very large x. Backwards compatibility has apparently resulted in these errors being perfectly preserved in all implementations of the ix87. Other ix87 errors are relatively minor AFAIK. E.g., exp(x) is not implemented directly in hardware. It requires several imprecise calculations like x*log2(e), so the hardware is good enough for double precision but not for long double precision where I guess it has an error of a few ulps. It has to switch the mode to long double precision to be good enough for double precision. log(x) is implemented more directly in hardware and I guess it has an error of about 1 ulp in long double precision. Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20080307163246.J22709>