Date: Tue, 1 Mar 2005 18:04:55 -0500 From: David Schultz <das@FreeBSD.ORG> To: Bruce Evans <bde@zeta.org.au> Cc: "M. Warner Losh" <imp@bsdimp.com> Subject: Re: cvs commit: src/lib/msun/src e_expf.c Message-ID: <20050301230455.GB3437@VARK.MIT.EDU> In-Reply-To: <20050302090138.K7588@delplex.bde.org> References: <200502240632.j1O6WDP9029589@repoman.freebsd.org> <20050226023149.GA63314@VARK.MIT.EDU> <20050226.180549.113100483.imp@bsdimp.com> <20050227023513.GA77813@VARK.MIT.EDU> <20050302090138.K7588@delplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
On Wed, Mar 02, 2005, Bruce Evans wrote: > On Sat, 26 Feb 2005, David Schultz wrote: > >>: This is related to a good reason why we can't switch the default > >>: precision on i386 to extended. Many of the functions in libm use > >>: minimax approximations, which are ``optimal'' approximations in > >>: the sense that their maximum error over all in-range inputs is the > >>: smallest possible (unless more terms are used). These approximations > >>: take rounding error into account, so when the machine precision is > >>: increased, they're no longer optimal and the error in the approximation > >>: can increase significantly. ... > > I didn't know that there was such a good reason. > > The library could switch to the appropriate mode. That wuld be > inefficient, but it may need to be done in some cases anyway when the > user has changed the FP environment. I'm not familiar with all the > details of what C99 allows. It allows or requires changes to the > rounding mode to affect most functions, but shouldn't this only affect > rounding of the result and not break algorithms that depending on > certain rounding internally? IIRC, changing the precision is much slower than changing the rounding mode. (Intel had to make changing the rounding mode fast because of C's broken float->int conversion semantics.) But yes, you're right that we could do this if we decide to change the rounding precision on x86. Ideally, libm functions should round the right way for the current rounding mode, although most of them probably don't do so. The polynomials they use internally, as well as the splitting tricks to get extra precision, probably only work for round to nearest. I was careful about setting the rounding mode appropriately in fma(), and it would not be hard to fix the rest of libm. It turns out that once you get round-to-nearest right, the other rounding modes are not very hard. This is because in round-to-nearest, you have to account for the possibility that the error term might be something like 0.50000000000001 * 2^-p, whereas the other rounding modes just require one extra bit of precision. In most of the libm routines, you'd probably just have to set the rounding to round-to-nearest, then change it back to whatever it was before immediately before the final calculation (often of the form a+b, where b is the tail of a). Very few people care about evaluating libm functions in those other rounding modes, so hopefully the performance impact of doing so would not be very large. > We support the extension of changing the > rounding precision. gcc does constant folding to 53-bit precision on FreeBSD/i386, so we clearly don't support it very well. > "small" (1-ulp) errors in these functions are 2^11 times larger than > 1-ulp errors in extended precision, so you get few benefits from > extended precision if you call one of these functions and it is > broken. Even if it isn't broken, then it will have 0.5-ulp errors. > Mainly code that never calls a math function or calls only long double > math functions that actually exist and actually are accurate to nearly > 1 extended precision ulp can benefit. Have you tested what the error is with extended precision? For some of the transcendental functions, I wouldn't be surprised if you get an error greater than 1 double precision ulp.
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20050301230455.GB3437>